Errors in the code when solving the transport equation

Question

I updated the content of the question, as I myself understood and corrected some points. And so there are the following variables:

hx (шаг по пространству)= 0.1 ht (шаг по времени) = 0.5 Nx (Количество шагов по пространству) = 10; Nt (Количество шагов по времени) = 12; wx[i] - массив содержит все шаги по пространству wt[j] - массив содержит все шаги по времени wht[j][i]-двумерный массив по которому будет строится результирующий график;

And so several problems were found:

1) Found in code

 for(int i = 0; i < Nx; i++) { wx[i+1] = wx[i] + hx; //массив wht[0][i] = fn(T, i * hx); //i * hx//Исправлено. }

where fn is:

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

I incorrectly set the initial conditions

Like wine on a chart, among the disgrace there is a common solution “wave” somewhere, but the rest is wrong there just because I incorrectly set the initial conditions in the code

2) Cycles and arrays, namely:

There is a call to a non-existent element of the array and therefore there are terrible numbers and a jump in the graph, but I cannot change the indices because there is a certain formula.

 for(int j = 0;j<Nt;j++) { for(int i = 1;i<Nx-1;i++)///Исправлено с помощью ответа пользователя Denis Protopopov(и график немного изменился) { 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]); } }

If I write instead of i = 0 , i = 2 or j = 2 , then the drawing of the chart will be terrible.

All code:

 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; } public:static double fn(int T,double x)//Функция Хэвисайда - Начальное условие { if (x>=0) return T; else if(x<0) return 0; } public: void drawfunc(double T) { double xmin = -5; double xmax = 10; for(double x = xmin;x<xmax;x+=0.01) { chart1->Series["Функция Хэвисайда"]->BorderWidth=3; chart1->Series["Функция Хэвисайда"]->Points->AddXY(x,fn(T,x)); } } private: System::Void button1_Click(System::Object^ sender, System::EventArgs^ e) { ///-------Переменные для работы с разностными схемами double a;//скорость переноса double hx,ht;//шаги по пространству и времени int Nx,Nt;//На сколько частей разбивае отрезок сетки int T;//Входной параметр для искомой функции(Хэвисайд) double wx[20]={0};//массив точек x double wt[20]={0};//массив точек t double wht[20][20]={0};// массив для разностной схемы(сетки) //double res[20][20]={0};//результирующий массив if(textBox1->Text!="" && textBox2->Text!="" && textBox3->Text!="" && textBox4->Text!="" && textBox5->Text!="" && textBox6->Text!="" && textBox6->Text!="") { ///-----Ввод данных----- a=Convert::ToDouble(textBox1->Text); hx=Convert::ToDouble(textBox2->Text); q=Convert::ToDouble(textBox3->Text); ht=Convert::ToDouble(textBox4->Text); Nx=Convert::ToInt32(textBox5->Text); Nt=Convert::ToInt32(textBox6->Text); T=Convert::ToInt32(textBox7->Text); //----Построение сетки,узлов 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]); } ///-----------------Вычисление основых формул с разностной схемой for(int j = 0;j<Nt;j++) { for(int i = 1;i<Nx-1;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])+ht*fn(wt[j],hx); } } //----Для записи в простой текстовый файл String^ fileName = "results.txt"; StreamWriter^ sw = gcnew StreamWriter(fileName); for(int j = 0;j<Nt;j++) { for(int i = 0 ;i<Nx;i++) { sw->Write("{0} ",wht[j][i]); } sw->WriteLine(); } sw->Close(); //Тестовые функции для отображения drawfunc(T);// Вызов функции для рисования ///---------------Построение графика for(int j = 0;j<Nt;j++) { for(int i = 0;i<Nx-1;i++) { chart2->Series["Series2"]->BorderWidth=3; chart2->Series["Series2"]->Points->AddXY(i,wht[j][i]); } } } else MessageBox::Show("Ошибка!Не все поля заполнены!"); }

At first glance, this is something related to mathematical analysis, or the geometry of graphs.
You already tell for the nesmyshlenny what exactly you are trying to achieve in this code?
The mathematical part is tested and the formula is purely programmed correctly (already checked) I want to achieve the correct output of the graph.
But problem number 1, prevents the conclusion of the schedule and therefore I am interested to know what actions I need to do to correctly set the initial conditions and work with the formula.
The fact is that I am trying to write the value of the fn function into an array where the steps over space are stored, which is fundamentally wrong.
The run must go from 0 to N-1 for (int i = 0; i <N-1; i ++) in order not to go beyond the array boundaries.
and tried to use the possibilities of commercial mat.paketov code generation?
First, the result of the first assignment does not affect the result, since the value of wht[j+1][i] overwritten by the second assignment.
Secondly, the second assignment is equivalent to the relation (wht[j+1][i] - wht[j][i])/ht == -a * ( wht[j][i+1] + wht[j][i]) / hx .
This, up to O (hx), is equivalent to hx * du/dt = -2au — a rigid differential equation.
Let me remind you that the transfer equation is du/dt = -a*du/dx .
Before you solve with a discontinuous initial condition, make sure that you have a smooth function - the sine, will run in the right direction at the right speed.

Accepted Answer · 2017-02-24T13:00:12

In order to answer why the wave does not run on the graph, you must first deal with the numerical method, and only then with the code in C ++. The theory of numerical analysis is not included in the range of issues that are discussed on this site, but brief information may be useful.

1. A small digression devoted to difference schemes

So, we are looking for the function u (x, t), which is a solution to the Cauchy problem for a partial differential equation

  dudu --- = - a --- (1) dtdx

with initial condition

  u(x,0) = f(x) (2)

Equation (1) is called the transport equation , since the solution to our problem is u(x,t) = f(x-at) . This is a function whose graph moves with time to the right at a speed of a .

Substitute the function f (x-at) instead of u in equation (1) and calculate the partial derivatives to verify that it really is a solution.

Teachers teach students not to look for easy ways, but to solve this equation using the finite difference method. We can offer several difference approximations of equation (1). We are invited to use one of them. As it turns out, this is a very bad choice.

So, consider the scheme in which the derivative with t to t approximated by the difference operator Dt u .

  Dt u = (u(x,t+dt) - u(x,t))/dt.

This operator has the first order of approximation.

We approximate the derivative with x to x difference operator Dx u

  Dx u = (u(x+dx,t) - u(x,t))/dx.

Will get

  Dt u = -a Dx u.

Specialists in numerical analysis are well aware that this scheme has the first order of approximation in dt and dx, and is also absolutely unstable . This means that any small disturbance this scheme increases exponentially.

We substitute the function q^(t/dt) e^(ipx/dx) instead of u(x,t) q^(t/dt) e^(ipx/dx) . Substituting, we see that if this function is a solution of a difference equation, then | q | > 1 if sin p/2 != 0 . That is, any small perturbation of such a form will grow as t grows. This method of studying the stability of difference schemes is called the von Neumann method .

So, this scheme is absolutely unstable, therefore it does not satisfy the condition of the Lax theorem, so one should not expect the convergence of the function calculated by this scheme to the solution of the initial problem for the partial differential equation.

Interestingly, if we calculate the derivative of x using the formula

  D'x u = (u(x,t) - u(x-dx,t))/dt,

then we get a stable scheme.

2. Statement of the problem on a limited interval

We want to search for a solution not on the whole numerical axis, but on a limited interval [xmin,xmax] . In order for the function f(x-at) be the solution to our problem, it is necessary to supplement the task with the following boundary conditions.

  u(xmin,t) = f(xmin - a*t) (3) u(xmax,t) = f(xmax - a*t) (4)

Now we can proceed to programming the method for solving the problem (1-4) based on explicit first-order approximation schemes.

3. Software implementation

Below is an example implementation of the method itself. Everything related to visualization is decided separately.

The first block of code is the definition of boundary conditions.

 #include <cstdio> #include <cstdlib> #include <cmath> using namespace std; const double T = 1.0; const double a = 0.5; inline double f(double x){ if (x>=0) return T; return 0; } inline double u0(double x){ return f(x); } inline double u1(double xmin, double t){ return f(xmin - a * t); } inline double u2(double xmax, double t){ return f(xmax - a * t); }

The following code block contains a function that calculates the solution and saves it to an array. Our function takes as its arguments the parameters of the problem - the coordinates of the boundaries of the calculated area and a link to the array in which to save the result.

 void run(double xmin,double xmax, double tmax, int N,double u[][M]){ double ht = tmax/N; double hx = (xmax - xmin)/M; for(int i=0; i<N;i++){ u[0][i] = u0(xmin + i*hx); } for(int j=0; j < N-1; j++){ u[j+1][0] = u1(xmin,j*ht); for(int i=1; i<N; i++){ double dudx = (u[j][i] - u[j][i-1])/hx; u[j+1][i] = -a * dudx * ht + u[j][i]; } } }

Here the variant with a stable scheme is implemented - the difference operator D'x

4. Results

In the following picture, the values of the function u(x,t) drawn for a=0.5 and t=0 , t=5.0 and t=10.0 .

It can be seen that the wave is not so much running as smeared in space.

Now we replace the main loop with a loop that implements an absolutely unstable circuit defined by the Dx operator.

  for(int j=0; j < N-1; j++){ for(int i=0; i<N-1; i++){ double dudx = (u[j][i+1] - u[j][i])/hx; u[j+1][i] = -a * dudx * ht + u[j][i]; } u[j+1][N-1] = u2(xmax,j*ht); }

The result will look something like this.

This result is not at all like a runner to the right.

5. Conclusion

I hope that this example illustrates that before looking for errors in the code, it is necessary to study the properties of the algorithm.

Looking at your code, I understand how I have complicated everything. And please tell me, and the 1st cycle shows instability.
xmax = as I understand it, Nx is my number where I will end the axis in space example 1 xmin = as I understand it 0 tmax = and I have this Nt (that is, axis maximum in time), for example, 1.2 N is the number Steps on the axes as I understand it?
I run run with xmin=-5.0 , xmax = 10.0 and tmax = Nt * ht .
The arrays wt and wh , are missing in my code - this is irrelevant here, and u replaces wht .
I almost guessed besides M. Already now it shows better, but not to the end.
The last question is the structure of the cycle for displaying a graph like this?
for (int j = 0; j <N-1; j ++) {for (int i = 1; i <N; i ++) {// Graph output}} // I use tChart.
I have one more cycle, the 3rd one, which is responsible for building a purely point of coordinates x.
In order not to fill the entire graph with an even layer, you should choose several values of j corresponding to different points in time, for example, j0 = 0, j1 = Nt / 2, j2 = Nt.
And for each of these moments, write the wht[ j0 ][i] values for all i=0..Nx to a separate TChartSeries instance.
The outer one runs through the list of selected indices j0, j1, j2.

GHosT GHosT 434 four 14 · Answer 2 · 2017-02-22T22:06:23

Since no one wants to deal with this issue, I dare to express their views.

The question is missing part of the code, namely:

No values are given to all elements of the wx , wt and wht arrays. When an array is declared, it is allocated a place in memory and, if the declaration does not specify values for the elements of the array, then when accessing them you get what it turned out to be. From here and these terrible values are taken. ))
Where is the variable T declared ??
Where is the variable a declared ??

Somewhat changed code:

 double hx = 0.1; // Шаг по пространству double ht = 0.5; // Шаг по времени int Nx = 10; // Количество шагов по пространству int Nt = 12; // Количество шагов по времени double wx[Nx]; // Массив содержит все шаги по пространству double wt[Nt]; // Массив содержит все шаги по времени double wht[Nt][Nx]; // Двумерный массив, по которому будет строится результирующий график for (int i = 0, step = 0; i < Nx; i++, step += hx) { wx[i] = step; // Заполняем массив. Каждый шаг больше предыдущего на hx. //wht[0][i] = fn(T, i * hx); — Непонятно назначение метода fn. Проще сделать так: wht[0][i] = (T < 0) ? 0 : i * hx; // В этом цикле заполняется только массив `wx` и первая строка массива `wht`, // а где задаются значения остальных элементов массивов `wht` и `wt` ?? } for (int j = 0; j < Nt; j++) { for (int i = 0; i < Nx; i++) { // Если в данном цикле использовать индекс j+1 или i+1, то это обязательно вызовет // ошибку, т.к. индекс выйдет за пределы массива. Ошибка как раз будет в конце // перебора массива. Думаю, этим и вызваны возникновение "страшных цифр" и скачок // графика. 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]); } }

@beginner In general, I would like to see the entire project as a whole, because of these pieces of code it is not possible to make a complete picture.