Программирую уравнение переноса

[beginer]

Pphantom
Спасибо за совет,я буду по ситуация применять оба ваших совета,зависит от постановки моей задачи

И всё-таки вопрос что мне делать дальше чтобы найти решение уравнения переноса?с правильным графиком.

выразил $y^{j+1}_{i}$а как дальше быть?Что нужно ещё найти или выразить мне,чтобы найти решение уравнения переноса и построить правильный график?

Дело в том, что я иду по шагово(для меня это оптимальный вариант), то есть нахожу одно,а потом второе а потом третье и так плавно к решению подхожу.

вот к примеру формула
$y^{j+1}_{i} = - \lambda \cdot y^{j}_{i+1}+(1+\lambda) \cdot y^{j}_{i}$
оказалась неправильной в моём случае.

Везде в книгах(У Самарского) есть метод гармоник получается что я должен его находить уже?

[Pphantom]

beginer в сообщении #1190921 писал(а):

И всё-таки вопрос что мне делать дальше чтобы найти решение уравнения переноса?с правильным графиком.

Формулу наконец-то написали правильную. Теперь осталось ее правильно запрограммировать, и получится то, что нужно (ну, с точностью до того, что явная схема первого порядка будет размазывать контактный разрыв, поэтому вместо движущейся "ступеньки" Вы будете получать хотя и перемещающийся нужной скоростью, но все более пологий перепад).

[beginer]

Я изменил код и пустил через цикл(по i и по j соответственно)но снова выдаёт совсем не то.
Сразу скажу что Nx=10,Nt=12(это мои количества отрезков на осях),ht=0,3 а hx=0,1

Вот сам код:

Используется синтаксис C++

for(int i = 0;i<Nx;i++)///Вывод
                                                 {
                                                         for(int j = 0;j<Nt;j++)
                                                         {
                                                                 wht[j+1][i]=(-a*ht*wht[j][i+1]+a*ht*wht[j][i]+wht[j][i]*hx)/hx;
                                                         }
                                                 }

НО я обратил внимание на то, что возможно я не так начальные\граничные условия задаю.
К примеру в Самарском(да и не только) начальные и граничные условия представлены в таких формулах.

$u(x,0)=u_0(x)$
$u(0,t) = u_0(t)$

Вот мой код

Используется синтаксис C++

 public:static double qx0(double x)//Начальные условия
                   {
                           if (x<=0) 
                                   return 0;
                           else
                                   return x;
                   }

        public:static double qt0(double t)//Граничные условия
                   {
                           if (t<=0) 
                                   return 0;
                           else
                                   return t;
                   }

//----Построение сетки,узлов

                                                 for(int i = 0;i<Nx;i++)        
                                                 {
                                                         wx[i+1]=wx[i]+hx;//массив 
                                                         wht[0][i]=qx0(wx[i]);
                                                 }      

                                                 for(int j = 0;j<=Nt;j++) //
                                                 {                                                       
                                                         wt[j+1]=wt[j]+ht;
                                                         wht[j+1][0] = qt0(wt[j+1]);                                                                                                                                                                             
                                                 }

А вот результат построения сетки(в текстовом виде,специально проверял)

И вот тут интересно,если было бы неправильно то такая сетка и не строилась(пусть даже в таком виде), но в тоже время есть такая мысль,что что-то неправильно сделано мной.

Поэтому интересно мнение знающих людей.



P.S
worm2 в этой теме написал прочитать характеристику уравнения переноса,я ознакомился - это общее решение этого уравнения. Но пока не применяю,так как и некуда.

[Pphantom]

beginer в сообщении #1191040 писал(а):

Я изменил код и пустил через цикл(по i и по j соответственно)но снова выдаёт совсем не то.

Естественно. У Вас снаружи идет цикл по координате, а внутри - по времени.

[beginer]

В результате следующего кода

Используется синтаксис C++

for(int j = 0;j<Nt;j++)
                                                 {
                                                         for(int i = 0;i<Nx;i++)///Вывод формул
                                                         {                                              
                                                                wht[j+1][i]=(-a*ht*wht[j][i+1]+a*ht*wht[j][i]+wht[j][i]*hx)/hx;
                                                         }
                                                 }

                                                for(int j = 0;j<Nt;j++)
                                                 {
                                                         for(int i = 0;i<Nx;i++)
                                                         {
                                                                 chart2->Series["Series2"]->BorderWidth=3;                      
                                                                 chart2->Series["Series2"]->Points->AddXY(i,wht[i][j]);
                                                         }
                                                 }      

Выводится такой график



Не похоже..и числа такие страшные,а задавал то скорость 4,шаги hx=0,1 ht=0,3,Nx=10,Nt=12.

Я уже и отдельно для координаты цикл добавил и всё равно не так.
Пока что смотрю в сторону начальных и граничных условий(где у меня 50 на 50)

[Pphantom]

beginer в сообщении #1191340 писал(а):

В результате следующего кода

... программа периодически читает что-то из wht[j][i+1] в тот момент, когда i+1=Nx. В общем, из несуществующего в массиве элемента. Что она оттуда извлечет, никому неведомо, соответственно, результат такого "моделирования" может быть каким угодно.

P.S. Кстати, использовать для вычислительного моделирования C# Вам кто-то подсказал? Или сами догадались... Просто я еще понимаю, когда это делается потому, что этот язык уже известен, а что-то другое учить лень/некогда/поздно. Но начинать с нуля... В общем, Вы какой-то уникальный образец человека, который умудряется делать все возможные ошибки (а иногда - и невозможные). Как Вам так удается?

[beginer]

Используемый язык программирования C++,который не изучается с нуля как и C#
Использую его потому что нужно по заданию.

А теперь по поводу кода.
Неизвестный элемент массива. Вся сетка инициализирована(см.рисунок с текстовым файлом),нули это тоже элементы массива и они существуют,несмотря на то что они 0.

Что касаемо цикла то я не спорю со людьми,которые умнее и знают больше меня,но я не могу понять тогда что же я делаю не так.

Объясняю почему не могу понять.

Вот есть разбиения на отрезки(количество шагов) Nt,Nx(12 и 10 соответственно).
Каждый шаг ht имеет свой слой.На каждом слое получается и будет вычисляться правильная формула
То есть иду вверх по оси t и на каждом шаге вычисляю каждый слой.

Вот что у меня происходит в коде(или что я представляю что у меня должно происходить)
-Цикл от j->после уже цикл по i выполняется 10 раз(решаем 10 раз правильную формулу),и только потом цикл от j прибавляется и так до конца(в данном случае до 12).
Но i отвечает то за ось x,следовательно поднялись по оси по времени на шаг вверх,и начинаем вычислять формулу при этом уже принимая по внимание ось x(шаги hx)

Вот что я ещё представляю.

$(ht_t,hx_{x1}),(ht_t,hx_{x2}),(ht_t,hx_{x3}),(ht_t,hx_{x4})$ и так далее..

При этом заметьте $ht_t$ этой как раз шаг по времени и он не изменяется,меняются только шаги по пространству.
То есть на каждом шаге по времени я выполняю 10 вычислений с координатами от шага по пространству.

[Pphantom]

beginer в сообщении #1191532 писал(а):

Используемый язык программирования C++,который не изучается с нуля как и C#

В общем-то в данном случае разница невелика.

beginer в сообщении #1191532 писал(а):

Использую его потому что нужно по заданию.

Печально.

beginer в сообщении #1191532 писал(а):

А теперь по поводу кода.
Неизвестный элемент массива. Вся сетка инициализирована(см.рисунок с текстовым файлом),нули это тоже элементы массива и они существуют,несмотря на то что они 0.

Увы, нет. Вы каждый раз на правом краю пространственной сетки используете значение на предыдущем временном слое, которое оказывается за границей массива. Происходит это потому, что для производной по координате везде используется правая аппроксимация, которую нельзя применять в самых правых точках сетки. Чтобы такое не происходило, вообще-то нужно принимать какие-то специальные меры. Либо для правого края сетки сменить аппроксимацию на левую, либо задать какие-то краевые условия (если их можно откуда-то взять).

Правда, в данной конкретной задаче лучше все-таки с самого начала поменять вычислительную схему. Использованная Вами разностная аппроксимация представляет собой т.н. "схему с разностью по потоку", которая, вообще говоря, неустойчива. Это, правда, ошибка уже "следующего порядка", но рано или поздно, когда Вы исправите более простые, Вы на нее наткнетесь. Равно как и на проблему с неусточивостью из-за слишком большого числа Куранта:

beginer в сообщении #1191340 писал(а):

Не похоже..и числа такие страшные,а задавал то скорость 4,шаги hx=0,1 ht=0,3

[beginer]

Текущее положение дел
Да с правой частью схемы буду работать(устойчивость,сходимость с помощью метода гармоники,уравнение пока решается и для этого и тему отдельную создал)

Попробовал схему с левой производной(как и тут посоветовали)
Формула:

$\frac{y^{j+1}_{i} - y^{j}_{i}}{\tau}-a \frac{y^{j}_{i} - y^{j}_{i-1}}{h}=0$

Её решение:
$y^{j+1}_{i}=\frac{a \tau y^{j}_{i} + a \tau y^{j}_{i-1} + y^{j}_{i} h}{h}$
Код:

Используется синтаксис C++

 for(int j = 0;j<Nt;j++)
                                         {
                                                 for(int i = 1;i<Nx;i++)/
                                                 {                      
                                                                
                                                         wht[j+1][i]=(-a*ht*wht[j][i]+a*ht*wht[j][i]+wht[j][i-1]*hx)/hx;                                                        
                                                 }
                                         }

И вышел вот такой график(конечно же он неправильный)



Как сказали на проверке проблема в начальных и граничных условиях,в качестве краевых условий решили испробовать единичную функцию Хэвисайда(попробовать в начальных условиях или граничных условиях(по очереди и не вместе конечно же))

Вот пытаюсь собственно с помощью функции Хэвисайда решить данную проблему
Код и график

Используется синтаксис C++

public:static double fn(double T,double x)//Функция Хэвисайда
                   {
                           if (x>=0)
                                   return T;
                           else if(x<0)
                                   return 0;  
                   }

И задаю её в в коде

Используется синтаксис C++

 for(int i = 0;i<Nx;i++)        
                                         {
                                                 wx[i+1]=wx[i]+hx;//массив 
                                                 wht[i][0]=qx0(wx[i]);
                                                // wht[i][0]=fn(i*hx,T);
                                         }      

                                         for(int j = 0;j<Nt;j++) //
                                         {                                                       
                                                 wt[j+1]=wt[j]+ht;
                                                 wht[0][j+1]= fn(T,wt[j]);
                                         }

И я думаю что в данном коде не задаю никаких граничных начальных условий вообще..

[Pphantom]

beginer в сообщении #1192600 писал(а):

И я думаю что в данном коде не задаю никаких граничных начальных условий вообще..

Вообще-то граничные условия - это одно, начальные условия - это другое. Что Вы задаете в коде, неведомо, поскольку полный код, кажется, ни разу не появлялся.

И, да, как Вы представляете себе граничные условия в виде функции Хевисайда? Не с точки зрения задания таковых (это так или иначе не проблема), а с точки зрения осмысленности подобного действия?

beginer в сообщении #1192600 писал(а):

уравнение пока решается и для этого и тему отдельную создал

Кстати, было бы лучше этого не делать. Конечно, кому как, но мне уже несколько надоело вспоминать, в какой теме мы обсуждали, как умножать $2\times 2$, а в какой - $2\times 3$.

[beginer]

Функция отвечающая за граничное условие(Единичная Хэвисайда)
$T$ сделал постоянной(временно),а сама функция с одним параметром $t$

Используется синтаксис C++

public:static double fn(double x)//Функция Хэвисайда - Начальное условие(Граничное условие) а начальное 0-1
                   {
                           double T=5;
                           if (x>=0)
                                   return T;
                           else if(x<0)
                                   return 0;  
                   }

График функции строится правильно.

А вот так отображаю это в цикле(задаю начальное условие)

Используется синтаксис C++

 for(int i = 0;i<Nx;i++)        
                                         {
                                                 wht[i][0]=fn(i*hx);
                                         }

а вот основной график на текущий момент



Ну и конечно цикл подсчёта и цикл графика вместе с ним

Используется синтаксис C++

for(int j = 0;j<=Nx;j++)
                                         {      
                                                 for(int i = 0;i<=Nt;i++)///Вывод формул wht[j+1][i] проследить за циклами
                                                 {                                                      
                                                         wht[i][j+1]=((-a*ht*wht[i][j]+a*ht*wht[i][j]+wht[i-1][j]*hx)/hx);
                                                 }
                                         }

                                         for(int j = 0;j<=Nt;j++)
                                         {
                                                 for(int i = 0;i<=Nx;i++)
                                                 {
                                                         chart2->Series["Series2"]->BorderWidth=3;                      
                                                         chart2->Series["Series2"]->Points->AddXY(i,wht[i][j]);
                                                 }
                                         }      

Но вот что я заметил

Я к примеру хочу чтобы координата шла по-шагово (и отображалось по шагово) а у меня как видно всё упирается в одну точку от неё идёт и обратно же в неё возвращается и следовательно такой ужасный график

[beginer]

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


Как мне сказали условие правильности будет тогда,когда верхний слой будет 0,следующий слой будет 1 и потом уже возрастают числа.

У меня же здесь одни пятёрки(Что конечно же неправильно)
Цифра 5,потому что я в Единичной функции Хэвисайда задаю константно 5 и потом уже пытаюсь в коде инициализировать начальные условия

[beginer]

Так,вроде одна из ошибок найдена.

Простая забывчивость приоритетов арифметических операций в C++ способно привести к страшным результатам.

[beginer]

Вот и плавно подошёл к 2-ум основным проблемам которые мне мешают построить мой график

1)Задание начальных условий .

Используется синтаксис C++

for(int i = 0;i<Nx;i++)        
                                         {
                                                 wx[i+1]=wx[i]+hx;//массив 
                                                 wht[0][i]=fn(T,i*hx);
                                         }  

Всё же какие действия нужно совершить мне чтобы задать их правильно?
Какие действия мне нужно совершить чтобы сделать код с правильным заданием граничных условий?

2) Неясность с индексами массивов в формуле

Да как отмечал уважаемый Pphantom в формуле идёт обращение к несуществующему элементу массива и следовательно возникают страшные цифры и скачок графика

Используется синтаксис C++

for(int j = 0;j<Nt;j++)
                                         {
                                                  for(int i = 0;i<Nx;i++)///Вывод формул wht[j+1][i] проследить за циклами
                                                 {                                                      
                                                         wht[j+1][i]=((wht[j+1][i]-wht[j][i])/ht)+a*((wht[j][i+1]-wht[j][i])/hx);
                                                         wht[j+1][i]=-a*(ht*(wht[j][i+1] + wht[j][i])/hx) + (wht[j][i]); //Линия идёт туда
                                                 }
                                         }

Но стоит мне например задать в цикле не i=0,j=0, а i=2,j=2 например то график становится совсем другим.

Пока что промежуточный график такой.



Как видно среди безобразия на графике,видно как прячется там подобие волны(решение уравнения переноса)

[beginer]

К предыдущему посту добавлю такую информацию что
wx[i]- содержит шаги по оси по пространству
Nx - это количество шагов по пространству(количество разбиений на отрезки)
hx - это значение шага(0,1)

а wht[j][i] это массив который и выводит результат

Взять к примеру этот код

Используется синтаксис C++

for(int i = 0;i<Nx;i++)        
                                         {
                                                 wx[i+1]=wx[i]+hx;//массив 
                                                 wht[0][i]=fn(T,i*hx);
                                         }  

Что я подразумеваю под работой этого кода?

Я записываю в массив шагов по пространству значение функции Хэвисайда и следовательно у меня в i записываются значения функции(учитывая что функция возвращает число по оси t,то это число и заполняет у меня всю ось пространства)

И это явно не может быть заданием начального условия.