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

[beginer]

Продолжение темы topic114556.html я и спрашивал про сетки и про двумерный массив пытался что-то сказать

1)Задаю исходные параметры шаги по пространству,шаг по времени,скорость переноса,концентрат-q,количество разбиений по x(Nx),количество разбиений по thau(Nt)

Нахожу начальные условия
$q(x,t)$\left\lvert$_{t=0} = q_{0}(x)$
и граничные условия(которые были даны мне для выполнения)
$q(x,t)$\left\lvert$_{x=0}=q_{0}(t)$
Код:

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

public:static double nach(double q,double hx)//Начальные условия
{
     return q*hx;
}

public:static double konc(double q,double ht)//Граничные условия
{
     return q*ht;
}
 

2) Как известно решение данного уравнения это $u_{0}(x-a\cdot t) $ .Следовательно добавляем отдельно функцию,учитывая уже известные начальные и граничные условия.

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

public:static double Psi(double hx,double ht,double a)///Функция решения уравнения переноса
                   {
               double qq;
               return nach(qq,hx)-(a*konc(qq,ht));
                   }

А теперь моменты в которых у меня возник вопрос.
Согласно моему заданию мне нужно использовать для начала
явную разностную схему

$\frac{y^{j+1}_{i} - y^j_{i}}{\tau}$+a$\frac{y^j_{i+1} - y^j_{i}}{h}=0$ (формула 1)

$y^{j+1}_{k}=$\gamma$ \cdot y^{j}_{k+1}+(1+$\gamma$) \cdot y^{j}_{k}$(формула 2) , где $\gamma=\frac{a\tau}{h}$


Мой примерный алгоритм:

1)Для начала я нахожу число Куранта

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

gama=(a*ht)/hx;//Число Куранта

2)После нахождения числа Куранта запрограммировал формулу-1 и формулу-2

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

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

                                         for(int i = 1;i<Nx+1;i++)      
                                         {
                                                 for(int j =1 ;j<Nt+1;j++)
                                                 {
                                                         wht[i][j]=gama*wht[j][i+1] + (1 + gama)*wht[j][i] + wht[j][i];//формула 2                                        
                                                 }      
                                         }

Пока значение wht у меня нули-нулями и поэтому вопрос
1)Правильный ли алгоритм?
2)Функция решения названное Psi(см. в коде) куда её применять?
Дело в том что $y^{j}_{i}$ обозначается и как сеточная функция и как отдельный массив где будут размещены полученные значения уравнений.

[beginer]

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

Поэтому буду спрашивать по частям
1)Начну с задания граничных и начальных условий

Вместо q во многих учебниках пишут u,но суть не важна.
Вот эти 2 формулы.
$q(x,t)$$\left\lvert_{t=0} = q_{0}(x)$
и граничные условия(которые были даны мне для выполнения)
$q(x,t)$$\left\lvert_{x=0}=q_{0}(t)$


Что же я сам подразумеваю по этими формулами - q нулевое зависит от $\tau$ которое у меня ht(шаг по времени) и зависит от x(hx - шаг по пространству и сетке).

Иными словами это $y=f(x)$ но в более другой формулировке и не забываем про стремящиеся к нулю $(x = 0,t = 0)$.
После повторения теории я решил для того чтобы найти начальные и граничные условия по указанным формулам нужно в $q_{0}$ подставить значения вышеперечисленных шагов и умножить их на $q$- концентрат.

В итоге я ввожу q-концентрат c клавиатуры(передаю в функции параметры в виде шагов и концентрата) и по отдельности получаю начальные и граничные условия.

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

public:static double nach(double q,double hx)//Начальные условия
{
     return q*hx;
}

public:static double konc(double q,double ht)//Граничные условия
{
     return q*ht;
}

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

[worm2]

1) Неправильно обозначать одной и той же буквой ($q_0$) и начальное, и граничное условия. Хотя бы потому, что это функции от разных переменных. Пусть будет $q_{x0}$ и $q_{t0}$. Думаю, что это основная причина неразберихи.

2) Зачем-то тут $\tau$, оно же $h_t$. Лучше пусть будет $h_t$ и $h_x$, а $\tau$ не нужно.

3) $q_{x0}$ зависит от $x$, а не от $h_x$. Так же $q_{t0}$ зависит от $t$, а не от $h_t$. $h_x$ и $h_t$ — это шаги по пространственной переменной $x$ и по времени $t$, в данной задаче это постоянные величины.

4) Ваши функции nach и konc не годятся. Я так понимаю, здесь имеются в виду граничные и начальные условия ($q_{x0}$ и $q_{t0}$ соответственно). Тогда у этих функций должен быть один аргумент. Например:

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

double qx0(double x)
{
  return sin(x); // Тут должно стоять граничное условие из вашего задания
}

double qt0(double t)
{
  return cos(t); // Тут должно стоять начальное условие из вашего задания
}

5) Вопрос на засыпку. Что не так в выдуманной мной паре условий ($q_{x0}(x)=\sin x$, $q_{t0}(t)=\cos t$) ?

[beginer]

Спасибо за подсказку а теперь я перехожу к следующему шагу-вопросу

Как известно решение данного уравнения это

$$u_{0}(x-a\cdot t) $$

Волновое.

Я решил отдельно вынести эту функцию в код,учитывая уже известные начальные и граничные условия.

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

public:static double funcsolve(double x,double t,double a)///Функция решения уравнения переноса
                   {
                           return qx0(x)-(a*qt0(t));
                   }

Почему я так сделал?
Обычно в программах которые настроены на решение диф.уравнений любого вида(обязательно встречается f,которая использует искомую функцию для дальнейших расчётов и решений по ней)

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

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

2)Возможно должно быть так
$$q_{x0}(x)=\sin t, $q_{t0}(t)=\cos x) $ так как при x t=0 $а при $t$ следовательно x = 0$

[worm2]

beginer писал(а):

$$u_{0}(x-a\cdot t) $$

Волновое.

Да, это общее решение уравнения переноса. Только я не помню, чтобы кто-то называл это решение волновым. Оно, действительно, похоже на волну тем, что распространяется с постоянной скоростью $a$. Но тут может возникнуть путаница с волновым уравнением. И ещё, я бы тут нолик убрал. Лишний он. $u$ — это и есть искомая функция, $u_0$ могла бы быть, например, функцией из граничных условий, а тут нолик сбивает с толку.

Итак, мы знаем общее решение уравнения переноса. Как его теперь получить из начального и граничного условий? Подставить сначала $x=0$, а потом $t=0$. Помня, что $q_{x0}$ задано только при $x\geqslant0$, а $q_{t0}$ задано только при $t\geqslant0$.

Цитата:

2)Возможно должно быть так $q_{x0}(x)=\sin t$, $q_{t0}(t)=\cos x$

Нетушки, если мы пишем (x), значит, зависит от $x$, а не от $t$. Вот попробуйте для моего примера найти решение $u$. Найдите в учебнике определение понятия "характеристика" для уравнения переноса.

-- Ср янв 18, 2017 20:01:26 --

Цитата:

1)Как я знаю в начальных\граничных условиях возможно экспоненциальное поведение...

Тоже неправильно. С точки зрения математики, возможно всякое, точнее не всякое, а математика тут накладывает свои ограничения, но они не такие.

[beginer]

И так используя функцию

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

public:static double funcsolve(double x,double t,double a)///Функция решения уравнения переноса
                   {
                           return x -(a*t);
                   }

получился такой график(уже почти похож на бегущую волну)


Оранжевая линия это функция,представленная кодом выше,а синяя это метод Рунге-Кутта(зачем он нужен это объясню позже)

А теперь перехожу такому моменту.
Согласно заданию даны следующие формулы
$\frac{y^{j+1}_{i} - y^j_{i}}{\tau}+a\frac{y^j_{i+1} - y^j_{i}}{h}=0$ (формула 1)

$y^{j+1}_{k}=$\gamma$ \cdot y^{j}_{k+1}+(1+$\gamma$) \cdot y^{j}_{k}$(формула 2) , где $\gamma=\frac{a\tau}{h}$ - это Число Куранта.

Для начала нахожу Число Куранта

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

gama=(a*ht)/hx;//Число Куранта

А теперь перехожу к формуле явной разностной схемы(формула 1)

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

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

И формула-2 это формула один после аппроксимации 1-ого порядка

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

for(int i = 1;i<Nx+1;i++)      
                                         {
                                                 for(int j =1 ;j<Nt+1;j++)
                                                 {
                                                         res[i]=gama*wht[j][i+1] + (1 + gama)*wht[j][i] + wht[j][i];//формула 2                                        
                                                 }      
                                         }

После нескольких опытов и полученных странных значений возник вопрос

Куда мне подставить решение моей функции в формулы-1 и 2(не код)?
Дело в том что $y^j_{i}$ может быть так и двумерным массивом или же функцией решения уравнения переноса.
Если сказать честно то вот здесь я запутался.

[beginer]

На основе вводимых параметров


и следующего кода

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

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++) //N+1 j+1
                                         {
                                                 wt[j+1]=wt[j]+ht;
                                                 wht[j+1][0] = qt0(wt[j+1]);
                                                 for(int i = 1;i<Nx+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]=gama*wht[j][i+1] + (1 + gama)*wht[j][i] + wht[j][i]+ht*fn(wx[i],wt[j]); 
                                         }
                                 }
 

Получаю сначала исходную матрицу(с отображением осей x t с шагами)

После применения 2 формул получаю такие значения


Потом я буду исследовать её на устойчивость и Рунге-Кутт(но это потом)
Вновь интересно мнение,правильно ли я всё делаю.

-- 22.01.2017, 23:42 --

Единственное что я сам заметил так это сдвиги(смотрим последний блокнотовский файл)

И следовательно каждая строка это возможно очередной слой

[beginer]

Сформулировалась такая задача по перенос

С помощью уравнения переноса нужно Единичную Функцию Хэвисайда "перенести"
Результат у меня должен быть сдвиг по оси и отрисовка функции Хэвисайда несколько раз(иными словами перенос по оси)

Функцию Хэвисайда запрограммировал вот так

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

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

Результат выполнения кода


Но дальше возник вопрос

Данные формулы нужно использовать
$\frac{y^{j+1}_{i} - y^j_{i}}{\tau}+a\frac{y^j_{i+1} - y^j_{i}}{h}=0$ (формула 1)

$y^{j+1}_{k}=\gamma \cdot y^{j}_{k+1}+(1+\gamma) \cdot y^{j}_{k}$(формула 2) , где $\gamma=\frac{a\tau}{h}$

В коде выглядит это так

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

 gama=(a*ht)/hx; //Число Куранта
               for(int i = 0;i<Nx;i++) 
                           {
                        for(int j = 0;j<Nx;ij++)
                                {
                                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]=gama*wht[j][i+1] + (1 + gama)*wht[j][i];
                                }
                        }

Разумеется результат на графике показывает что никаких сдвигов и переноса нету.

1)В каком направлении нужно мне двигаться чтобы прийти к правильному ответу?
2)Функцию Хэвисайда нужно в любом случае кроме вывода на экран куда-то применять и мне её подставлять в формулу где разностная схема или нужно выразить какой-либо коэффициент?

[Pphantom]

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

$y^{j+1}_{k}=\gamma \cdot y^{j}_{k+1}+(1+\gamma) \cdot y^{j}_{k}$(формула 2) , где $\gamma=\frac{a\tau}{h}$

Не надо использовать эту формулу, она неверна.

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

В коде выглядит это так

Простите, но это чушь. Даже если отвлечься от попыток реализации неправильной модели.

Я Вам это уже писал, но повторюсь, и учтите, что это не преувеличение, а суровая реальность: Вам нужно пока забыть про "уравнения переноса" и тому подобные сложные вещи и изучить нормально школьную математику 7-9 классов. А заодно, как выясняется, и программирование (школьное или первого курса, где оно у Вас там было). Пытаться получить осмысленный результат, случайным образом комбинируя действия, смысл которых Вы не понимаете и не пытаетесь понять, невозможно.

[beginer]

Спасибо за ответ про формулу.
Значит будут другие попытки и выложу сюда

P.S(ОФФТОП)
А вот по поводу вспомнить,я вспоминаю медленно но верно,но мне главное это сдать,а как это сдам уже вспомню многое и без нервов

[beginer]

Пока вот что получилось



Для правильности решения нужно выровнять и убрать скачки которые остаются внизу

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

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

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

где fn это моя функция Хэвисайда,а wx[i]-моя ось с шагами h=0.1, а wt[j]-ось t с шагами \tau

Пока мысль такова что я просто не те параметры передаю в свою функцию fn..

[beginer]

В результате работы следующего кода,график функции немного изменился.

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

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

Как оказалось в функции fn нужно было другие параметры,а именно оригинальные параметры fn(T,x). А я просто вставил оси x,t.
Да график немного стал получше..
Но до "ступенек" ему ещё далеко

[Pphantom]

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

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

Она опять неверна.

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

В результате работы следующего кода,график функции немного изменился.

Вот у Вас две пары вложенных циклов. Вы можете объяснить - зачем?

[beginer]

Это было в качестве эксперимента(как выяснилось 2 цикла по отдельности делать это неудачное решение)

На данный момент я выразил $y^{j+1}_i$ (Здесь на форуме мне указали мои ошибки,спасибо)
Вот полная формула(Она была проверена руководителем и верна,процесс решения уравнения записан аналитически)

$y^{j+1}_i=\frac{-a\cdot \tau y^j_{i+1} + a\cdot \tau y^j_{i} + y^j_{i+1} \cdot h}{h}$

Вопрос
до нахождения $y^{j+1}_i$

Получается я должен найти аппроксимацию,сходимость,устойчивость формулы явной разностной схемы и только потом выразить $y^{j+1}_i$?

Я прошу простой ответ наподобие да или нет?

P.S.
В учебниках Самарского и не только,в большинстве случаев идёт сначала формула схемы,а потом и последующие действия такие как аппроксимация,нахождения сходимости и устойчивости
но ключевое слово в большинстве(так как бывает сразу в метод гармоник уходит вычисление)

[Pphantom]

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

Это было в качестве эксперимента(как выяснилось 2 цикла по отдельности делать это неудачное решение)

Делать их вместе - решение не менее неудачное.

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

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

$y^{j+1}_i=\frac{-a\cdot \tau y^j_{i+1} + a\cdot \tau y^j_{i} + y^j_{i+1} \cdot h}{h}$

Ну слава богам. В этот раз действительно правильно.

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

Получается я должен найти аппроксимацию,сходимость,устойчивость формулы явной разностной схемы и только потом выразить $y^{j+1}_i$?

Я прошу простой ответ наподобие да или нет?

Ха. Из этих двух вариантов, по-видимому, ближе "нет", но в целом вопрос почти лишен смысла.

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