Numerical solution of system of PDEs (1 order,linear,homogenious) method Lax-Wendroff

Question

I'm trying to solve linear system of PDEs $$\frac{\partial \bar{u}}{\partial t}+\begin{pmatrix} \frac{2}{5} & -\frac{72}{5} & \frac{64}{5} \\ \frac{6}{5} & \frac{-56}{5} & \frac{42}{5} \\ \frac{12}{5} & \frac{-52}{5} & \frac{34}{5} \\ \end{pmatrix} \frac{\partial \bar{u}}{\partial x}=0$$ with conditions

• $u_{1}(x,0)=3,x<0$
• $u_{1}(x,0)=1,x>0$
• $u_{2}(x,0)=1,x<0$
• $u_{2}(x,0)=2,x>0$
• $u_{3}(x,0)=2,x<0$
• $u_{3}(x,0)=3,x>0$

numerically (with a help of method of Lax-Wendroff). So, I've written a program on MATLAB.

function [Xg,Yg,st,APlus,AMinus] = Sett21(N)  
  A=[2/5 -72/5 64/5 ; 6/5 -56/5 42/5 ; 12/5 -52/5 34/5]  
  A1=zeros(3,3);  
  for i=1:3  
    for j=1:3  
      if (A(i,j)>0)  
        A1(i,j)=A(i,j);  
      else  
        A1(i,j)=0;  
      end;  
    end;  
  end;  
  A2=zeros(3,3);  
  for i=1:3  
    for j=1:3  
      if (A(i,j)<0)  
        A2(i,j)=A(i,j);  
      else  
        A2(i,j)=0;  
      end;  
    end;  
  end;  
  for i=1:N  
    x(i)=-3+(6/N)*(i-1);  
  end;  
  for i=1:(10*N)  
    y(i)=(5/(10*N))*(i-1);  
  end;  
  Xg=x;  
  Yg=y;  
  st=3;  
  APlus=A1;  
  AMinus=A2;

That makes a grid on (X,T) and something more, but this won't be used this time.

Then, I've made another program that calculates the values of $\bar{u}(x,t)$ in the points of the grid.

[Xg,Yg,st,APlus,AMinus] = Sett21(50);
N=50;
Ukt=zeros(3,N,10*N);
for t=1:3
  for i=1:N
    if (t==1)
      if (Xg(i)<=0)
        Ukt(1,i,1)=3;
      else
        Ukt(1,i,1)=1;
      end;
    elseif (t==2)
      if (Xg(i)<=0)
        Ukt(2,i,1)=1;
      else
        Ukt(2,i,1)=2;
      end;
    elseif (t==3)
      if (Xg(i)<=0)
        Ukt(3,i,1)=2;
      else
        Ukt(3,i,1)=3;
      end;
    end;
  end;
end;
h1=6/N;
h2=5/(10*N);
A=[2/5 -72/5 64/5; 6/5 -56/5 42/5; 12/5 -52/5 34/5];
T1=[-2/5+h1/h2 72/5 -64/5; -6/5 56/5+h1/h2 -42/5; -12/5 52/5 -34/5+h1/h2];
T2=[2/5+h1/h2 -72/5 64/5; 6/5 -56/5+h1/h2 42/5; 12/5 -52/5 34/5+h1/h2];
for j=2:10*N
  for i=2:N-1
    vec=(1/2)*(eye(3)+(h2/h1)*A)*[Ukt(1,i-1,j-1);Ukt(2,i-1,j-1);Ukt(3,i-1,j-1)]+(1/2)*(eye(3)-(h2/h1)*A)*[Ukt(1,i+1,j-1);Ukt(2,i+1,j-1);Ukt(3,i+1,j-1)];
    Ukt(1,i,j)=vec(1);
    Ukt(2,i,j)=vec(2);
    Ukt(3,i,j)=vec(3);
  end;
  vec1=-A*[Ukt(1,2,j); Ukt(2,2,j); Ukt(3,2,j)]+(h1/h2).*[Ukt(1,1,j-1); Ukt(2,1,j-1);Ukt(3,1,j-1)];
  vecx=T1\vec1;
  Ukt(1,1,j)=vecx(1);
  Ukt(2,1,j)=vecx(2);
  Ukt(3,1,j)=vecx(3);
  vec2=A*[Ukt(1,N-1,j);Ukt(2,N-1,j);Ukt(3,N-1,j)]+(h1/h2).*[Ukt(1,N,j-1);Ukt(2,N,j-1);Ukt(3,N,j-1)];
  vecxx=T2\vec2;
  Ukt(1,N,j)=vecxx(1);
  Ukt(2,N,j)=vecxx(2);
  Ukt(3,N,j)=vecxx(3);
end;
for i=1:N
  for j=1:10*N
    Ukt1(i,j)=Ukt(1,i,j);
  end;
end;
for i=1:N
  for j=1:10*N
    Ukt2(i,j)=Ukt(2,i,j);
  end;
end;
for i=1:N
  for j=1:10*N
    Ukt3(i,j)=Ukt(3,i,j);
  end;
end;
mesh(Yg,Xg,Ukt1);

Something goes not right and the solution between characteristics of this system reaches infinity. I know that it is not right, because it is the Reimann's problem, that's why the solution is constant between characteristics. I've been searching for a mistake since previous month and can't find it. Please, help me.

Answer 1

N=400;

Ukt=zeros(3,N,10*N);
for i=1:N
x(i)=-3+(6/N)*(i-1);
end;
for i=1:(10*N)
 y(i)=(5/(10*N))*(i-1);
 end;

Xg=x;
Yg=y;
for t=1:3
   for i=1:N
       if (t==1)
          if (Xg(i)<=0) Ukt(1,i,1)=3;
          else Ukt(1,i,1)=1;
            end;
       elseif (t==2)
             if (Xg(i)<=0) Ukt(2,i,1)=1;
             else Ukt(2,i,1)=2;
             end;
       elseif (t==3)
            if (Xg(i)<=0) Ukt(3,i,1)=2;
            else Ukt(3,i,1)=3;
            end;
       end;
    end;
end;

Answer 2

1st mistake - APlus and AMinus are defined incorrectly. $A^{+}=R\Lambda^{+}L; A^{-}=R\Lambda^{-}L$, where $\Lambda^{+},\Lambda^{-}$ are positive and negative parts of matrix with eigenvalues, $R,L$ are matrix with right and left eigenvectors accordingly.

2nd mistake - it is not Lax-Wendroff method at all. That has absolutely another formula.

Answer 3

h1=6/N;
h2=5/(10*N);
alfa=h2/h1;
A=[2/5 -72/5 64/5; 6/5 -56/5 42/5; 12/5 -52/5 34/5];
T1=[-2/5+h1/h2 72/5 -64/5; -6/5 56/5+h1/h2 -42/5; -12/5 52/5 -34/5+h1/h2];
T2=[2/5+h1/h2 -72/5 64/5; 6/5 -56/5+h1/h2 42/5; 12/5 -52/5 34/5+h1/h2];
for j=1:10*N-1
for i=2:N-1
vec=(eye(3)-(alfa^2)*A^2)*[Ukt(1,i,j); Ukt(2,i,j); Ukt(3,i,j)]+((1/2)*alfa*A+(1/2)*   (alfa^2)*A^2)*[Ukt(1,i-1,j); Ukt(2,i-1,j); Ukt(3,i-1,j)]+((-1/2)*alfa*A+(1/2)*(alfa^2)*A^2)*[Ukt(1,i+1,j);Ukt(2,i+1,j);Ukt(3,i+1,j)];
Ukt(1,i,j+1)=vec(1);
Ukt(2,i,j+1)=vec(2);
Ukt(3,i,j+1)=vec(3);
end;
vec1=-A*[Ukt(1,2,j+1); Ukt(2,2,j+1); Ukt(3,2,j+1)]+(h1/h2)*[Ukt(1,1,j);   Ukt(2,1,j);Ukt(3,1,j)];
vecx=T1\vec1;
Ukt(1,1,j+1)=vecx(1);
Ukt(2,1,j+1)=vecx(2);
Ukt(3,1,j+1)=vecx(3);
vec2=A*[Ukt(1,N-1,j+1);Ukt(2,N-1,j+1);Ukt(3,N-1,j+1)]+(h1/h2)*[Ukt(1,N,j);Ukt(2,N,j);Ukt(3,N,j)];
vecxx=T2\vec2;
Ukt(1,N,j+1)=vecxx(1);
Ukt(2,N,j+1)=vecxx(2);
Ukt(3,N,j+1)=vecxx(3);
end;
 for i=1:N
 for j=1:10*N
    Ukt1(i,j)=Ukt(1,i,j);
 end;
 end;
 for i=1:N
 for j=1:10*N
    Ukt2(i,j)=Ukt(2,i,j);
 end;
 end;
   for i=1:N
   for j=1:10*N
    Ukt3(i,j)=Ukt(3,i,j);
 end;
 end;
    %mesh(Yg,Xg,Ukt3)