2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


В этом разделе нельзя создавать новые темы.

Если Вы хотите задать новый вопрос, то не дописывайте его в существующую тему, а создайте новую в корневом разделе "Помогите решить/разобраться (М)".

Если Вы зададите новый вопрос в существующей теме, то в случае нарушения оформления или других правил форума Ваше сообщение и все ответы на него могут быть удалены без предупреждения.

Не ищите на этом форуме халяву, правила запрещают участникам публиковать готовые решения стандартных учебных задач. Автор вопроса обязан привести свои попытки решения и указать конкретные затруднения.

Обязательно просмотрите тему Правила данного раздела, иначе Ваша тема может быть удалена или перемещена в Карантин, а Вы так и не узнаете, почему.



Начать новую тему Ответить на тему
 
 Вычислительная математика
Сообщение24.05.2013, 15:24 
Аватара пользователя


12/05/12
604
Оттуда
Здравствуйте. Помогите, пожалуйста, найти ошибку. Это - программа на языке 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}$- шаги по координате и по времени. Но увы, программа выдаёт разрушенное решение. Разрушение происходит именно между разрывов, там, где, по идее, должно быть постоянное решение. То есть, ошибка в исходной формуле. Уже месяц не могу найти эту ошибку. Товарищи,помогите, пожалуйста.

 Профиль  
                  
 
 Re: Вычислительная математика
Сообщение24.05.2013, 16:18 
Заслуженный участник


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

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

 Профиль  
                  
 
 Re: Вычислительная математика
Сообщение24.05.2013, 16:31 
Аватара пользователя


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

 Профиль  
                  
 
 Re: Вычислительная математика
Сообщение25.05.2013, 02:29 
Заслуженный участник


23/07/08
10626
Crna Gora
Код:
vecx=T1\vec1;
Что это за операция?

Код:
(h1/h2).*[вектор, заданный компонентами]
Что означает такое умножение с точкой?

 Профиль  
                  
 
 Re: Вычислительная математика
Сообщение25.05.2013, 05:21 
Заслуженный участник
Аватара пользователя


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

Цитата:
Там реализована формула $\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}$

 Профиль  
                  
 
 Re: Вычислительная математика
Сообщение25.05.2013, 11:21 
Аватара пользователя


12/05/12
604
Оттуда
svv
Здесь вычислялись неявно крайние узлы.

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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 6 ] 

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group