偏微分各差分格式MATLAB編程

1、（1）迎風(fēng)格式function u = peHypbYF(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);if a>0 for j=1:(n+M) u0(j) = IniU(minx+(j-M-1)*h); endelse for j=1:(n+M) u0(j) = IniU(minx+(j-1)*h); endend u1 = u0; for k=1:M if a>0 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

2、 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.*(x<=0.5)+0.*(x>0.5);plot(x,y,'g')title('Ó­·ç¸ñʽ')format short;運行結(jié)果>> u=peHypbYF(1

3、,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 = u1;end x=0.01:0.01:n/100; u = u1

4、(M+1):(M+n);plot(x,u,'r'); axis(0 1 -0.5 1.5);hold ony=1.*(x<=0.5)+0.*(x>0.5);plot(x,y,'g')title('Lax-Fridrichs¸ñʽ')format short;運行結(jié)果>> u = LaxFridrichs(1,0.005,101,0,1,100)（3）Lax-Wendroff格式function u = LaxWendroff(a,dt,n,minx,maxx,M)format lo

5、ng;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)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

6、.*(x<=0.5)+0.*(x>0.5);plot(x,y,'g')title('Lax-Wendroff¸ñʽ')format short運行結(jié)果>> u= 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.*(x<=0.5)+0.*(x>0.5);plot(x,y,'g')title