0
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – user642796
    May 24, 2013 at 17:47
  • $\begingroup$ 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. $\endgroup$
    – user642796
    May 24, 2013 at 17:47
  • $\begingroup$ 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. $\endgroup$
    – cool
    May 24, 2013 at 18:32
  • $\begingroup$ 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.) $\endgroup$
    – user642796
    May 24, 2013 at 19:12
  • $\begingroup$ 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! $\endgroup$
    – cool
    Jun 2, 2013 at 17:24

3 Answers 3

0
$\begingroup$

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.

$\endgroup$
0
$\begingroup$
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;
$\endgroup$
0
$\begingroup$
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)
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .