Vectorizing to speed up integration

Hi everyone --

In the code below I am performing some symbolic manipulation of a function that I subsequently want to integrate. I am trying to stay symbolic as long as I can so that I do not have to repeat any calculations. Once I get close to my final goal, I want to pass the function an array of variables and do numerical intagration many times. This is where I am stuck. Through searching the forum and the MATLAB help I have seen a few suggestions involving either anonymous functions or further symbolic manipulation, but I just can't get it to work. I suspect someone with more experience could do it rather quickly. FWIW, I will need to perform these calculations on scores (or hundreds) of times on data sets that are >10k elements long. Thanks in advance.

% This is intended to be a function where rx, ry, rz, eps1, and eps2 are 
% passed as arguments. Each variable is Nx1 so the function returns Nx2
% values. I have hard-coded some variables here for illustration.

clear;

lo = 0.8;
hi = 1.2;
o = ones(4,1);

rx      = (lo+(hi-lo)*rand(4,1)).*o;
ry      = (lo+(hi-lo)*rand(4,1)).*o;
rz      = (lo+(hi-lo)*rand(4,1)).*o;
eps1    = (lo+(hi-lo)*rand(4,1)).*o;
eps2    = (lo+(hi-lo)*rand(4,1)).*o;

% This would then be the beginning of the function.

syms a b c e1 e2 real positive
syms u v real

% Parametric equation of a superellipsoid
r = [a.*(cos(u).^e1).*(cos(v).^e2), b.*(cos(u).^e1).*(sin(v).^e2), c.*sin(u).^e1];

% Partial derivatives
r_u = diff(r, u);
r_uu = diff(r_u, u);
r_v = diff(r, v);
r_vv = diff(r_v, v);
r_uv = diff(r_u, v);

% Normal vector to the surface
n = cross(r_u, r_v);
n = n / norm(n);

% First fundamental form
E = dot(r_u, r_u);
F = dot(r_u, r_v);
G = dot(r_v, r_v);

% Second fundamental form
L = dot(r_uu, n);
M = dot(r_uv, n);
N = dot(r_vv, n);

% Area element
A = E*G - F^2;
dA = sqrt(A);

% Gaussian curvature
K = (L*N - M^2) / A;

% Mean curvature
H = (L*G - 2*M*F + N*E) / (2*A);

% Curvature elements to integrate
dK = K * dA;
dH = H * dA;

% Some options here, none of which I could get to work very well:
% DK = matlabFunction(dK)          % syntax problems
% dk = int(int(dK,u,0,u),v,0,v)     % super slow
% @(...) integral2(@(...) ...,a,b,c,d)
% and so on....

% REFS
% https://www.mathworks.com/matlabcentral/answers/1978659-how-to-calculate-the-mean-integrated-curvature-of-an-ellipsoid
% https://www.mathworks.com/matlabcentral/answers/2172786-average-curvature-of-a-closed-3d-surface

Answer 1

[Walter Roberson]

You can improve speed a little by using

temp = int(int(dK,u,0,u,hold=true),v,0,v,hold=true)
dk = release(temp);

However, using matlabFunction() or integral2() or the like cannot work, as your expression is over the seven variables [a, b, c, e1, e2, u, v] but numeric approaches cannot handle unassigned variables (other than the variables of integration.) Also, you are doing indefinite integration (upper bounds are symbolic) rather than definite integration (upper bounds numeric.)

Symbolic integration is your only option.

[Moe Szyslak]

Thanks so much for the response. I will try this approach. When I tried just using int(int(...)), my slow laptop chugged along for nearly an hour before I aborted, though. That said....

In looking at my code, I can see that I was a bit (a lot?) hasty (read: sloppy) in what I provided for the final task. Ultimately, I want to integrate those two curvature functions over the surface originally defined by 'r'. Since it is symmetric and centered on the origin, it suffices to integrate for 0<=u,v<=pi/2 and multiply the answer by 8. I would like to do this many, many times for different values of (a,b,c,e1,e2).

[Walter Roberson]

F = matlabFunction(dK, 'vars', [u, v, a, b, c, e1, e2], 'file', 'DK.m');

%then later
G = @(u, v) F(v, u, A, B, C, E1, E2);       %specific numeric A B C E1 E2
result = integral2(G, 0, pi/2, 0, pi/2)

[Moe Szyslak]

Thank you again for the response. I tried something similar to this prior to posting the original question, but I did not save the function as a file, I just used it in place. However, the approach that you've suggested seems to have the same problem as that which I tried. Specifically, I would like to be able to pass the final function Nx1 arrays of A B C E1 E2 and have the integration be performed N times, but I always get an array dimension mismatch error when I try this. I suppose that the naive solution is to wrap the integration in a for loop, which is what I was trying to avoid to begin with.

Thanks so much for the suggestions.

[Torsten]

Since the integrations for different values (A,B,C,E1,E2) are independent, there is no other way than to use a for-loop (or the "arrayfun" function which is a for-loop in disguise).

[Moe Szyslak]

Thank you -- I appreciate the explanation.

[Walter Roberson]

integral2() internally passes 2D arrays of the independent variables to the integrand. In order for N x 1 arrays of A B C E1 E2 to have a chance of working, they would have to be 1 x 1 x N arrays. But unfortunately that simply isn't supported.

[Moe Szyslak]

Good to know -- thank you for the clarification.

[Moe Szyslak]

OK, sorry, I just noticed something in Walter's post using matlabFunction():

Is switching the order of u, v intentional (if so, why?) or a typo? Sorry, I am really just learning to jump between symbolics and anonymous functions.

[Torsten]

Since integrations with respect to u and with respect to v are both from 0 to pi/2, the order of u and v doesn't matter. But maybe it's less confusing if the order is retained.

F = @(u,v)3*u.^2+v;
G = @(u,v)F(v,u);
result = integral2(G,0,pi/2,0,pi/2)

result = 8.0260

G = @(u,v)F(u,v);
result = integral2(G,0,pi/2,0,pi/2)

result = 8.0260

[Walter Roberson]

it was a typo. I had switched them at one point because I had misread the constraints as if u was bounded by v, but a further reading showed that instead u had a lower bound of zero and is unbounded on the upper range. I then reread and realized that you were probably intending to bound both of them. anyway I switched back the order but must have missed one of the places.