格式试验已知初值当x0时u1当x0u0用三

1、uu0 ，初值，当 x 0 时， u1，当 x 0,u=0,用三种格式。格式试验：已知utx解：（ 1），采用 laxwendroff 格式：题设中方程常常写为守恒型：ux( u 2)0（ 1）t2这样可以得到如下方程：uF (u)0(2)tx方程 (2)的 laxwendroff差分格式为：u nj1u njt (F jn 1F jn 1 )t 22 Anj 1/ 2 ( Fjn 1 F jn ) AJn1/ 2 ( F jnF jn 1 )2x2 x对于 F=u*u/2,A=u,可令Ajn1/ 21/ 2(unj1u nj )Anj 1/ 21/ 2(u njunj 1 )以下用 Matl

2、ab 编程实现上述思想：程序如下，其中后面为解释语句。clear alldx=1;dt=0.6*dx;ve(1:20)=1;%Laxwendoff 格式 dt/dx=0.6 时ve(21:40)=0;plot(ve,'-bo')%Lax 格式解析解ve1(1:20)=1;%Lax 格式 dt/dx=1 时ve1(21:40)=0;%for i=1:timestepfor j=2:39ve(j)=ve(j)-1/4*0.6*(ve(j+1)+ve(j-1)*(ve(j+1)-ve(j-1)+1/8*0.36*(ve(j+1)+ve(j)2*(ve(j+1)-ve(j)-(ve(j

3、)+ve(j-1)2*(ve(j)-ve(j-1);ve1(j)=ve1(j)-1/4*(ve1(j+1)+ve1(j-1)*(ve1(j+1)-ve1(j-1)+1/8*(ve1(j+1)+ve1(j)2*(ve1(j+1)-ve1(j)-(ve1(j)+ve1(j-1)2*(ve1(j)-ve1(j-1);end%endhold onplot(ve,'-b.')%Lax 格式hold on;plot(ve1,'-rs')text(0,0.2,'Laxwendoff 格式的数值解');legend('EXACT SOLUTION'

4、;,'dt/dx=0.6','dt/dx=1',-1)数值模拟结果：从结果格式结果可以看到，这一格式色散强于耗散。（ 2），采用 lax 格式：同上，方程 (2)的 lax 差分格式为：u nj 11/ 2(u nj 1u nj 1 )t ( F jn 1 F jn 1 )2 x程序实现过程如下：%Lax格式%clcclear alltimestep=1;dx=1;dt=0.6*dx;velocity(1:20)=1;%Lax格式 dt/dx=0.6 时velocity(21:40)=0;plot(velocity,'-bo')%Lax格式解析解

5、velocity1(1:20)=1;%Lax格式 dt/dx=1 时velocity1(21:40)=0;%for i=1:timestepfor j=2:39velocity(j)=1/2*(velocity(j+1)+velocity(j-1)-dt/dx*(velocity(j+1)+velocity(j-1)*(velocity(j+1)-velocity(j-1)/4;velocity1(j)=1/2*(velocity1(j+1)+velocity1(j-1)-(velocity1(j+1)+velocity1(j-1)*(velocity1(j+1)-velocity1(j-1)

6、/4;end%endhold onplot(velocity,'-b.')%Lax格式hold on;plot(velocity1,'-rs')text(0,0.2,'Lax 格式的数值解');legend('EXACT SOLUTION','dt/dx=0.6','dt/dx=1',-1)数值模拟结果：从模拟结果中我们可以看出， 其耗散项强于色散项。 因此在计算有强间断问题时， 容易通过间断。（ 3）采用 MacCormack 格式这是一种预估修正格式。对线性方程为u nj1u njc(u nj

7、1 u nj )unj11 / 2(u jnu nj1c(u jn 1unj 11 )在预估步使用前差，在修正步则使用后差。也可先后差然后前差。这对某些问题有优越性。对于线性方程，它等价于L-W 格式。对于于（2）方程可写为：u nj1unjt ( F jn1 F jn )xu nj11/ 2u nju nj1t ( F jn 1F jn11 )x程序实现过程如下：clear alldx=1;dt=0.6*dx;ve(1:20)=1;%MacCormack格式 dt/dx=0.6 时ve(21:40)=0;plot(ve,'-bo')%Lax格式解析解ve1(1:20)=1;%

8、Lax格式 dt/dx=1 时ve1(21:40)=0;%for i=1:timestepfor j=2:39tempve1=ve(j)-1/2*0.6*(ve(j+1)+ve(j)*(ve(j+1)-ve(j); tempve2=ve(j-1)-1/2*0.6*(ve(j)+ve(j-1)*(ve(j)-ve(j-1); ve(j)=1/2*(ve(j)+tempve1-1/2*0.6*(tempve1+tempve2)*(tempve1-tempve2); tempve1=ve1(j)-1/2*(ve1(j+1)+ve1(j)*(ve1(j+1)-ve1(j); tempve2=ve1(j-1)-1/2*(ve1(j)+ve(j-1)*(ve1(j)-ve1(j-1); ve1(j)=1/2*(ve1(j)+tempve1-1/2*(tempve1+tempve2)*(tempve1-tempve2);end%endhold onplot(ve,'-b.')%Lax格式hold on;plot(ve1