二阶TVD格式算出来答不到二阶，求教

wenjuncai1

按照附件图片中的方法想用对流方程测试下其阶数，u_t+u_x=0, 周期边界，初值取 sin(pi*x), matlab 程序如下, 但是跑一会数值解的波峰和波谷就会变平掉，精度也的达不到2阶，请教为什么会有这个结果，程序里面什么地方有问题呢？ 万分感谢

1. clc, clear
   
2. f = inline('u');
   
3. N = 200; dx = 2/N; x = 0:dx:2-dx; x = x';
   
4. dt = 0.5*dx; tend = 50; T = tend/dt;
   
5. u = sin(pi*x); u1 = zeros(N,1);
   
6. 
   
7. for j = 1:T
   
8.     t = j*dt;
   
9.     um1 = circshift(u,1); um2 = circshift(um1,1); up1 = circshift(u,-1);
   
10.     ubar = (f(u)-f(um1))./(u-um1); lam = dt/dx*ubar;
    
11.    
    
12.     Id1 = find(ubar>0); Id2 = find(ubar<=0);
    
13.    
    
14.     if ~isempty(Id1)
    
15.         r1 = (um1(Id1)-um2(Id1))./(u(Id1)-um1(Id1));
    
16.         phi1 = superbee(r1); Phi1 = circshift(phi1,-1);
    
17.         
    
18.         u1(Id1) = u(Id1)-dt/dx*(f(u(Id1))-f(um1(Id1)))-...
    
19.             1/2*lam(Id1).*(1-lam(Id1)).*(Phi1.*(up1(Id1)-u(Id1))-phi1.*(u(Id1)-um1(Id1)));
    
20.     end
    
21.    
    
22.     if ~isempty(Id2)
    
23.         r2 = (up1(Id2)-u(Id2))./(u(Id2)-um1(Id2));
    
24.         phi2 = superbee(r2); Phi2 = circshift(phi2,-1);
    
25.         
    
26.         u1(Id2) = u(Id2)-dt/dx*(f(up1(Id2))-f(u(Id2)))-...
    
27.             1/2*lam(Id2).*(1+lam(Id2)).*(Phi2.*(up1(Id2)-u(Id2))-phi2.*(u(Id2)-um1(Id2)));
    
28.     end
    
29.    
    
30.     uex = sin(pi*(x-t));
    
31.     plot(x,u1,x,uex,'r--');  h = gca; set(h,'xlim',[0,2],'Ylim',[-1.5,1.5]);
    
32.     drawnow
    
33.    
    
34.     u = u1;
    
35.     fprintf('Computing times = %d      Total times = %d\n',j,T)
    
36.    
    
37. end
    
38. 
    

复制代码