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

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.

-
It has been suggested that this question be migrated to one of our sister sites on the SE network, namely Computational Science. In my opinion, it is still within math.SE's domain, although I could not say where you are more likely to get a better answer. If you wish this to be migrated there, please respond to this comment and say so. – arjafi May 24 '13 at 17:47
Regardless of where this question ends up you may wish to improve its formatting as it is quite difficult to make out. Perhaps peruse our Markdown Help page for some further information on these matters. – arjafi May 24 '13 at 17:47
Perhaps it would be better if this question had been migrated to Computational Science, because, as I see, there are not so many mathematicians who know numerical methods. – cool May 24 '13 at 18:32
I have attempted to better format the question. Please double-check to make sure no mistakes have been introduced. Once you are happy with this, I'll migrate it to Computational Science. (Note that you will also want to make an account on Computational Science, if you haven't already.) – arjafi May 24 '13 at 19:12
Thanks, I have already found 2 mistakes here, first - in defining matrices APlus and AMinus (they must be defined with a help of eigenvalues) and in formula of Lax-Wendroff. Now everything works! – cool Jun 2 '13 at 17:24

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.

-
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;

-
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)

-