ep2_p5.m

此代码展示了Euler计算方法与差分方法

%  Ep2_p5: << 实验二 >>  参数方程的曲面图形表面积
%  Designed by FGH

m= 260;
R= 30;
[v,u]= meshgrid(0:pi/2/m:pi/2,0:2*pi/m:2*pi);
x= R*sin(v).*cos(u).*(1+0.1*abs(sin(6*u)));
y= R*sin(v).*sin(u).*(1+0.1*abs(sin(6*u)));
z= R*cos(v);

surfl(x,y,z,[-80,20]);
r= 0:1/256:1;
g= 0:1/256:1;
b= 0:0.3/256:0.3;
colormap([r' g' b']);
shading flat

clear all

syms u v R
%R= 30;

x= R*sin(v)*cos(u)*(1+0.1*sin(6*u));
y= R*sin(v)*sin(u)*(1+0.1*sin(6*u));
z= R*cos(v);

E= simple( diff(x,u)^2 + diff(y,u)^2 + diff(z,u)^2 );
G= simple( diff(x,v)^2 + diff(y,v)^2 + diff(z,v)^2 );
F= simple( diff(x,u)*diff(x,v) + diff(y,u)*diff(y,v) + diff(z,u)*diff(z,v) );
EG_F2= simple(E*G - F^2)
%f_uv= sqrt(EG_F2)
%f_u= int(f_uv,v,0,pi/2)
%f= int(f_u,u,0,pi/6)
%S= 12*double( f )

clear all

m= 15;
R= 30;
k= pi/6/(2*m);
h= pi/2/(2*m);
u= 0 : k : pi/6;
v= (0 : h : pi/2)';

 for j= 1:2*m+1
    for i= 1:2*m+1
       f(i,j)= sqrt( 1/100*R^2*sin(v(j))^2*(101+20*sin(6*u(i))+35*cos(6*u(i))^2)...
         *(1/100*R^2*cos(v(j))^2+1/5*R^2*cos(v(j))^2*sin(6*u(i))...
         -1/100*R^2*cos(v(j))^2*cos(6*u(i))^2+R^2)...
         -9/2500*sin(v(j))^2*R^4*cos(v(j))^2*cos(6*u(i))^2*(10+sin(6*u(i)))^2 );
   end
 end
 
clear Iij;
for j= 2:2:2*m
   for i= 2:2:2*m
%     算式(2.10)
      Iij(i,j)= k*h/9*( f(i-1,j-1)+f(i+1,j-1)+f(i-1,j+1)+f(i+1,j+1)...
              + 4*(f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1))...
              + 16*f(i,j) );
   end
end
I= sum(sum(Iij));
S= 12*I;
ans= sprintf(' 表面积%7.2f (m2), 需金箔%7.2f (m2)\n',S,(1+0.016)*S)
text(-50,-70,1,ans)


pause
close