clear;clc;close all
c=1000;
g=9.81;
lambda=0.025;
delta=0.005;
I=10; %x方向分段数
dx=0.2; %x方向的步长
T=0.2; %计算总时间
N=c; %t方向的分段数
dt=dx/c; %t方向的步长
A=1000;
w=1000;
%定义计算矩阵h和v
h=zeros(N+1,I+1);
v=zeros(N+1,I+1);
%计算h和v能确定的两个边界
%初始条件
h(1,:)=A;
v(1,:)=0;
t=0:dt:T;
%边界条件
h(:,1)=A*cos(w*t)';
v(:,I+1)=0;
NN=0;
eh=1; %定义误差
ev=1; %定义误差
msgbox('MATLAB编程答疑，请加QQ: 1530497909','MATLAB答疑','help')
web http://url.cn/TKcdXk -browser
while((ev>1e-5)|(eh>1e-5)) %当精度不满足时，继续迭代
h0=h; %保存前一次计算的h计算结果
v0=v; %保存前一次计算的v计算结果
for n=2:N+1 %计算h和v需要用内部点表示的边界
v(n,1)=(g/c)*(h(n,1)-h(n-1,2))+v(n-1,2)*(1-abs(v(n-1,2))*lambda*dt/(4*delta));
h(n,I+1)=h(n-1,I)+(c/g)*v(n-1,I)-v(n-1,I)*abs(v(n-1,I))*lambda*dx/(4*g*delta);
end
%计算中间点
for i=2:I
for n=2:N+1
h(n,i)=0.5*((h(n-1,i-1)+h(n-1,i+1))+(c/g)*(v(n-1,i-1)-v(n-1,i+1))-(lambda*dx/(4*g*delta))*(v(n-1,i-1)*abs(v(n-1,i-1))-v(n-1,i+1)*abs(v(n-1,i+1))));
v(n,i)=0.5*((g/c)*(h(n-1,i-1)-h(n-1,i+1))+(v(n-1,i-1)+v(n-1,i+1))-(lambda*dt/(4*delta))*(v(n-1,i-1)*abs(v(n-1,i-1))+v(n-1,i+1)*abs(v(n-1,i+1))));
end
end
%计算前后两个矩阵对应元素的最大误差
eh=max(max(abs((h0-h)./(h+eps))));
ev=max(max(abs((v0-v)./(v+eps))));
NN=NN+1; %记录循环次数
end
figure
subplot(3,1,1)
plot(t,h(:,1))
title('i=0处的h值随时间变化曲线')
subplot(3,1,2)
plot(t,h(:,6))
title('i=5处的h值随时间变化曲线')
subplot(3,1,3)
plot(t,h(:,11))
title('i=10处的h值随时间变化曲线')
[xx,tt]=meshgrid(0:I,t); %绘制x和y方向的网格线，以便绘制三维图
figure
mesh(xx,tt,h) %绘制h
xlabel('x')
ylabel('t')
zlabel('h')
figure
mesh(xx,tt,v) %绘制v
xlabel('x')
ylabel('t')
zlabel('v')
