P1_波浪对海上建筑物的作用_札记_色散方程解法（matlab）

《波浪对海上建筑物的作用》，电子版见评论区

两个色散方程，

 式2.18可用二分法或者matlab中一元函数零点的数值解法函数（fzero）

二分法示例（知k求w）

clc,clear,close all
k=0.1/20;
d=20;
g=9.81;
w=dispersion_w(k,d);
fprintf('k=%.5f,w=%.5f\n', k,w);

function w=dispersion_w(k,d)
w1=0;w3=100;i=0;g=9.81;  %设定w的上下限，一般w∈(0,5)
f=@(w)g*k*tanh(k*d)-w^2;
while (w3-w1>0.00001)  %设定进度
    w2mid=(w3+w1)/2;
    y1=f(w1);
    y2=f(w2mid);
    y3=f(w3);
    if y1*y2>0
        w1=w2mid;
    else
        w3=w2mid;
    end
    i=i+1;
end
w=w1;
end

式2.18、式2.20的 fzero 解法

%给定波数k0，求相应的w,再求出非传播模态对应的kn
k0=0.05; 
d=1; %水深
g=9.81;
f=@(w)g*k0*tanh(k0*d)-w^2;  %句柄函数
w=fzero(f,[0,100]);  %数值求零点的函数，fzero（函数，点附近或者某个区间内）

%solve evanescent kn
n=5;  %求出前n个解
kn=zeros(n,1);
f=@(k)g*k.*tan(k*d)+w^2;      %具有周期性的函数
for i=1:n
    a=(i-0.9999/2)*pi/d;   %设成1会有奇点
    b=(i+0.9999/2)*pi/d;
    kn(i,1)=fzero(f,[a,b]);
end

式2.20的图像特点，式2.20右侧图像如下，具有周期性，具有无穷多个解kn，且lim（kn+1-kn）=π/d，取水深d=π，图像代码见下方。

一般给定波数k，确定波频w，再由式2.20确定kn，式2.20右侧的图像是不会变动的，零点只会随着波数w变动，下午可直观看清左右两侧的值对应图像。

clc,clear,close all
d=pi;
k=0:0.005:10;
y1=-9.81*k.*tan(k*d);
w=3;
y2=w^2*ones(1,size(k,w));
plot(k,y1,'ob',k,y2,'om')
legend('y1','y2','FontSize', 24)
xlabel('kn')
ylim([-25,25])
grid on

书中的图像图2.1应有误，tan(kd)图像的周期性对不上，比如途中kd/π=4时，tan(4π)=0，但还是能够说明多根的，并且注意此处的k不是波数，而是一个普通变量。其下方的w^2*d/g=2是给定了一个工况，代码中给定是w=3，d=pi，有利于看清式2.20的周期性。