Vectorizing to speed up integration

1 次查看（过去 30 天）
Moe Szyslak
Moe Szyslak 2025-2-13
评论： Walter Roberson 2025-2-15
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

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2025-2-13
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.
  10 个评论
显示 8更早的评论
Torsten
Torsten 2025-2-15
编辑：Torsten 2025-2-15
Ran in:
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
Walter Roberson 2025-2-15

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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by