Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
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. –  Arthur Fischer 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. –  Arthur Fischer 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.) –  Arthur Fischer 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

3 Answers 3

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.

share|improve this answer
up vote 0 down vote accepted
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;
share|improve this answer
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)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.