michell积分程序.m

牛顿迭代法

format long
syms x z w
L=10;B=1.0;T=0.625;g=9.81;p=1000;        %L、B、T为主尺度；g为重力加速度；p为密度；分别采用国际单位
fx=eval(diff('(B/2)*(1-(2*x/L)^2)*(1-(z/T)^2)','x'));   %船体表面方程f对x的偏导
fz=eval(diff('(B/2)*(1-(2*x/L)^2)*(1-(z/T)^2)','z'));   %船体表面方程f对z的偏导
S1=inline(sqrt(1+fx^2+fz^2));            %积分S的被积函数表达式
S=2*dblquad(S1,-L/2,L/2,-T,0);           %船体静浮状态下的湿表面积S
i=0;
%计算不同航速下的Fn、Cw
for c=1.5:0.02:6;                        %c为船速 1.5m/s～6m/s
K=g/(c^2);                               %K为波数
i=i+1;
G1=(1-(2*x/L)^2)*(1-(z/T)^2)*(exp(K*(sec(w))^2*z))*cos(K*sec(w)*x);   %积分G的被积函数表达式
G2=int(G1,x,-L/2,L/2);                   %计算G的内层积分x
G3=int(G2,z,-T,0);                       %计算G的外层积分z 
G=B/2*G3;                                %求G
R1=inline((G^2)*(sec(w))^5);             %积分R的被积函数表达式
R=(4*p*g*(K^3)/pi)*quadl(R1,0,pi/2);     %求R
Cw(i)=0.5*R/(p*c^2*S);                   %不同航速下的兴波阻力系数
Fn(i)=c/sqrt(g*L);                       %不同航速下的弗氏数
end
plot(Fn,Cw);                             %绘出Fn-Cw曲线
xlabel('Fn'),ylabel('Cw');
title('兴波阻力系数随航速的变化');