Вычислительная математика

[cool.phenon]

Здравствуйте. Помогите, пожалуйста, найти ошибку. Это - программа на языке MATLAB для решения системы линейных однородных уравнений c частными производными 1 порядка методом Лакса-Вендроффа.

Создание сетки, нахождение положительной и отрицательной части матриц - остались из предыдущей программы.Здесь ошибок нет.

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;

А это само решение

[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=(eye(3)-(h2/h1)*(h2/h1)*A*A)*[Ukt(1,i,j-1); Ukt(2,i,j-1); Ukt(3,i,j-1)]+((h2/h1)*A+(1/2)*(h2/h1)*(h2/h1)*A*A)*[Ukt(1,i-1,j-1); Ukt(2,i-1,j-1); Ukt(3,i-1,j-1)]+((-h2/h1)*A+(1/2)*(h2/h1)*(h2/h1)*A*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);

Знаками восклицания выделена вероятная строка с ошибкой. Там реализована формула $\bar{u}_{k}^{n+1}=(E-\alpha^{2}A^{2})\bar{u}_{k}^{n}+(\alpha A+\frac{1}{2} \alpha^2 A^{2})\bar{u}_{k-1}^{n}+(-\alpha A+\frac{1}{2} \alpha^2 A^{2})\bar{u}_{k+1}^{n}$ - схема Лакса-Вендроффа для системы, $\alpha=\frac{h_{2}}{h_{1}};h_{2},h_{1}$- шаги по координате и по времени. Но увы, программа выдаёт разрушенное решение. Разрушение происходит именно между разрывов, там, где, по идее, должно быть постоянное решение. То есть, ошибка в исходной формуле. Уже месяц не могу найти эту ошибку. Товарищи,помогите, пожалуйста.

[svv]

Где ни читаю про эту схему, везде в итерационных формулах фигурирует $\frac{\Delta t}{\Delta x}$. А у Вас $\frac{h_2}{h_1}$, и Вы говорите, что $h_2$ и $h_1$ -- это шаги по координате и времени. То есть у Вас вроде как наоборот?

У Вас в программе местами $\frac{h_2}{h_1}$, местами $\frac{h_1}{h_2}$. Проверьте, во всех таких случаях записано то, что нужно?

[cool.phenon]

svv
Здесь просто неточность, конечно же, $h_{2}$ - шаг по времени, $h_{1}$ - шаг по координате. Те места, где используется $\frac{h_{1}}{h_{2}}$ - это вырезки неявных схем, там находятся самые крайние узлы на новых временных слоях. Там ошибки точно нету, потому что программу прогонял много раз, на краях условие должно быть постоянно, оно таким и остаётся.

[svv]

vecx=T1\vec1;

Что это за операция?

(h1/h2).*[вектор, заданный компонентами]

Что означает такое умножение с точкой?

[TOTAL]

Сформулируйте здесь задачу, которую решаете численно.
Что называете положительной (отрицательной) матрицей?

Цитата:

Там реализована формула $\bar{u}_{k}^{n+1}=(E-\alpha^{2}A^{2})\bar{u}_{k}^{n}+(\alpha A+\frac{1}{2} \alpha^2 A^{2})\bar{u}_{k-1}^{n}+(-\alpha A+\frac{1}{2} \alpha^2 A^{2})\bar{u}_{k+1}^{n}$

Проверьте эту формулу, надо

$\bar{u}_{k}^{n+1}=(E-\alpha^{2}A^{2})\bar{u}_{k}^{n}+(\frac{1}{2}\alpha A+\frac{1}{2} \alpha^2 A^{2})\bar{u}_{k-1}^{n}+(-\frac{1}{2}\alpha A+\frac{1}{2} \alpha^2 A^{2})\bar{u}_{k+1}^{n}$

[cool.phenon]

svv
Здесь вычислялись неявно крайние узлы.

TOTAL
Честное слово, Вы не представляете, как Вы меня выручили. Даже больше, чем выручили. Теперь у меня есть чем апеллировать к моему сенсею, потому что всё заработало. Но ошибка не моя, таков конспект.