How to calculate the mean integrated curvature of an ellipsoid?

7 ビュー (過去 30 日間)
Johan Zúñiga
Johan Zúñiga 2023 年 6 月 5 日
コメント済み: Johan Zúñiga 2024 年 9 月 14 日
I am using the following code to compute the surface area, mean integrated curvature, and Euler characteristic of the ellipsoid parameterized by r(x,y). Here 'x' and 'y' represent the angles $\theta$ and $\phi$ in spherical coordinates. The code is
clc; clear;
%ellipsoid radii
a=value1; b=value2; c=value2;
syms x y
%parametric vector equation of the ellipsoid
r=[a*cos(y)*sin(x), b*sin(x)*sin(y), c*cos(x)];
%partial derivatives
r_x=diff(r,x);
r_xx=diff(r_x,x);
r_y=diff(r,y);
r_yy=diff(r_y,y);
r_xy=diff(r_x,y);
%normal vector to the surface at each point
n=cross(r_x,r_y)/norm(cross(r_x,r_y));
%coefficients of the second fundamental form
E=dot(r_x,r_x);
F=dot(r_x,r_y);
G=dot(r_y,r_y);
L=dot(r_xx,n);
M=dot(r_xy,n);
N=dot(r_yy,n);
A=E*G-F^2;
K=((E*N)+(G*L)-(2*F*M))/A;
H=(L*N-M^2)/A;
dA=sqrt(A);
dC=K*dA;
dX=H*dA;
%converting syms functions to matlab functions
dA_=matlabFunction(dA); dC_=matlabFunction(dC); dX_=matlabFunction(dX);
% numerical integration
S = integral2(dA_, 0, 2*pi, 0, pi); %total surface area
C = (1/2)*integral2(dC_, 0, 2*pi, 0, pi); % medium integrated curvature
X = (1/(2*pi))*integral2(dX_, 0, 2*pi, 0, pi); % Euler characteristics
S and X give me very well, but the problem is with C, which always gives me very small, of the order of 1e-16 to 1e-14, which is absurd since C must be of the order of 4*pi*r (for the case of a sphere a=b=c=r). Could someone tell me what is wrong?

採用された回答

Milan Bansal
Milan Bansal 2024 年 9 月 11 日
編集済み: Milan Bansal 2024 年 9 月 13 日
Hi Johan Zúñiga
To calculate the mean integrated curvature of an ellipsoid, you need to evaluate the integral of the mean curvature HHH over the surface of the ellipsoid. For an ellipsoid, the mean curvature H at any point is the average of the two principal curvatures and
For a general surface, the mean curvature H is given by:
In the case of an ellipsoid, the mean integrated curvature can be computed by:
For a sphere with radius rrr, the mean integrated curvature is well-known to be:
For an ellipsoid with semi-axes a, b, and c, the exact analytical solution is not as straightforward as for the sphere. However, you can compute it numerically.
  • Mean Curvature H: Using the elements of the first and second fundamental forms, calculate the mean curvature H at each point on the ellipsoid surface using:
  • Numerically Integrate: Integrate the mean curvature H over the surface of the ellipsoid using integral2.
Here is how you can modify your code to calculate the curvature correctly:
clc; clear;
% Ellipsoid radii
a = 2; b = 5; c = 5;
syms x y
% Parametric vector equation of the ellipsoid
r = [a*cos(x)*sin(y), b*sin(x)*sin(y), c*cos(y)];
% Partial derivatives with respect to theta and phi
r_x = diff(r, x);
r_xx = diff(r_x, x);
r_y = diff(r, y);
r_yy = diff(r_y, y);
r_xy = diff(r_x, y);
% Normal vector to the surface
n = cross(r_x, r_y) / norm(cross(r_x, r_y));
% First fundamental form coefficients
E = dot(r_x, r_x);
F = dot(r_x, r_y);
G = dot(r_y, r_y);
% Second fundamental form coefficients
L = dot(r_xx, n);
M = dot(r_xy, n);
N = dot(r_yy, n);
% Determinant of the first fundamental form (area element factor)
A = E*G - F^2;
K=((E*N) + (G*L) - (2*F*M)) / A;
X=(L*N - M^2) / A;
% Mean curvature
H = (L*G - 2*M*F + N*E) / (2*A);
% Area element (sqrt(A))
dA = sqrt(A);
dC = K*dA;
dX = X*dA;
% Mean curvature element to integrate
dH = H * dA;
%converting syms functions to matlab functions
dA_ = matlabFunction(dA);
dC_ = matlabFunction(dC);
dX_ = matlabFunction(dX);
dH_ = matlabFunction(dH, 'Vars', [x, y]);
% Numerical integration over the parameter domain (theta = 0 to 2*pi, phi = 0 to pi)
S = integral2(dA_, 0, 2*pi, 0, pi) %total surface area
S = 200.0445
C_H = integral2(dH_, 0, 2*pi, 0, pi)
C_H = 52.3037
X = (1/(2*pi))*integral2(dX_, 0, 2*pi, 0, pi) % Euler characteristics
X = 2.0000
Hope this helps!
  1 件のコメント
Johan Zúñiga
Johan Zúñiga 2024 年 9 月 14 日
Thank you very much. Excellent corrections and adjustments.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeSurface and Mesh Plots についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by