MATLAB实现偏微分方程的差分计算 源程序代码.zip

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')