format long
syms x z w
L=10;B=1.0;T=0.625;g=9.81;p=1000;        %L、B、T為主尺度；g為重力加速度；p為密度；分別采用國際單位
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);     %求R
Cw(i)=0.5*R/(p*c^2);                   %不同航速下的興波阻力系數
Fn(i)=c/sqrt(g*L);                       %不同航速下的弗氏數
end
plot(Fn,Cw);                             %繪出Fn-Cw曲線
xlabel('Fn'),ylabel('Cw');
title('興波阻力系數隨航速的變化');