偏微分各差分格式MATLAB编程.doc

《偏微分各差分格式MATLAB编程.doc》由会员分享，可在线阅读，更多相关《偏微分各差分格式MATLAB编程.doc（22页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date偏微分各差分格式MATLAB编程偏微分各差分格式MATLAB编程（1）迎风格式function u = peHypbYF(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);if a0 for j=1:(n+M) u0(j) = IniU(minx+(j-M-1)*h); endelse for j=1:(n+M)

2、 u0(j) = IniU(minx+(j-1)*h); endend u1 = u0; for k=1:M if a0 for i=(k+1):n+M u1(i) = -(dt/h)*a*(u0(i)-u0(i-1)+u0(i); end else for i=1:n+M-k u1(i) = -(dt/h)*a*(u0(i+1)-u0(i)+u0(i); end end u0 = u1;endx=0.01:0.01:n/100; u = u1(M+1):(M+n);plot(x,u,r);axis(0 1 -0.5 1.5);hold ony=1.*(x0.5);plot(x,y,g)tit

3、le()format short;运行结果 u=peHypbYF(1,0.005,101,0,1,100)（2）Lax-Fridrichs格式function u = LaxFridrichs(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M) u0(j) = IniU(minx+(j-M-1)*h);end u1 = u0; for k=1:M for i=k+1:n+2*M-k u1(i) = -(dt/h)*a*(u0(i+1)-u0(i-1)/2+(u0(i+1)+u0(i-1)/2; end u0

4、= u1;end x=0.01:0.01:n/100; u = u1(M+1):(M+n);plot(x,u,r); axis(0 1 -0.5 1.5);hold ony=1.*(x0.5);plot(x,y,g)title(Lax-Fridrichs)format short;运行结果 u = LaxFridrichs(1,0.005,101,0,1,100)（3）Lax-Wendroff格式function u = LaxWendroff(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M) u0(j)

5、= IniU(minx+(j-M-1)*h);end u1 = u0; for k=1:M for i=k+1:n+2*M-k u1(i) = (dt/h)2*a*a*(u0(i+1)-2*u0(i)+u0(i-1)/2- . (dt/h)*a*(u0(i+1)-u0(i-1)/2+u0(i); end u0 = u1;endx=0.01:0.01:n/100; u = u1(M+1):(M+n);plot(x,u,r); axis(0 1 -0.5 1.5);hold ony=1.*(x0.5);plot(x,y,g)title(Lax-Wendroff)format short运行结果 u

6、= LaxWendroff(1,0.005,101,0,1,100)（4）Beam-Warming格式function u =BeamWarming(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M) u0(j) = IniU(minx+(j-M-1)*h);end u1 = u0; for k=1:M for i=k+2:n+2*M-k u1(i) = -a*(dt/h)*(u0(i)-u0(i-1)-a*(dt/h)*(1-a*(dt/h)*(u0(i)-2*u0(i-1)+u0(i-2)/2+u0(i); end u0 = u1;end x=0.01:0.01:n/100; u = u1(M+1):(M+n);plot(x,u,r); axis(0 1 -0.5 1.5);hold ony=1.*(x0.5);plot(x,y,g)title(Beam-Warming)format short;运行结果 u= BeamWarming(1,0.005,101,0,1,100)-