Ошибки в коде при решении уравнения переноса

Question

Обновил содержание вопроса,так как некоторые моменты я сам понял и исправил. И так есть следующие переменные:

hx (шаг по пространству)= 0.1

ht (шаг по времени) = 0.5

Nx (Количество шагов по пространству) = 10;

Nt (Количество шагов по времени) = 12;

wx[i] - массив содержит все шаги по пространству

wt[j] - массив содержит все шаги по времени

wht[j][i]-двумерный массив по которому будет строится результирующий график;

И так найдены несколько проблем:

1) Найдена в коде

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

где fn это:

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

Я неправильно задаю начальные условия

Как вино на графике, среди безобразия там где-то видно общее решение "волна", но остальное там неправильно как раз из за того что я в коде неправильно задаю начальные условия

2) Циклы и массивы, а именно:

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

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]);
    }
}

Если я напишу вместо i = 0, i = 2 или j = 2, то отрисовка графика будет ужасной.

Весь код:

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("Ошибка!Не все поля заполнены!");
             }

Answer 1

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

1. Небольшое отступление, посвященное разностным схемам

Итак, мы ищем функцию u(x,t), которая является решением задачи Коши для дифференциального уравнения в частных производных

  d u       d u
  --- = - a ---       (1)
  d t       d x

с начальным условием

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

Уравнение (1) называется уравнением переноса, поскольку решением нашей задачи является u(x,t) = f(x-at). Это функция, график которой движется со временем вправо со скоростью a.

Подставьте функцию f(x-at) вместо u в уравнение (1) и вычислите частные производные, чтобы проверить, что она действительно является решением.

Преподаватели учат студентов не искать лёгких путей, а решить это уравнение методом конечных разностей. Можно предложить несколько разностных аппроксимаций уравнения (1). Нам предлагается воспользоваться одной из них. Как вскоре выяснится, это очень неудачный выбор.

Итак, рассмотрим схему, в которой производная по t аппроксимируется разностным оператором Dt u.

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

Этот оператор обладает первым порядком аппроксимации.

Производную по x аппроксимируем разностным оператором Dx u

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

Получим

    Dt u = -a Dx u.

Специалистам по численному анализу хорошо известно, что эта схема обладает первым порядком аппроксимации по dt и dx, а также является абсолютно неустойчивой. Это означает, что любое малое возмущение эта схема увеличивает в геометрической прогрессии.

Подставим вместо u(x,t) функцию q^(t/dt) e^(ipx/dx). Подставив, увидим, что если эта функция является решением разностного уравнения, то |q| > 1, если sin p/2 != 0. То есть, любое малое возмущение такой формы будет расти при росте t. Этот способ исследования устойчивости разностных схем называется методом фон Неймана.

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

Интересно, что если мы будем вычислять производную по x по формуле

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

то получим устойчивую схему.

2. Постановка задачи на ограниченном интервале

Мы хотим искать решение не на всей числовой оси, а на ограниченном интервале [xmin,xmax]. Для того, чтобы функция f(x-at) была решением нашей задачи, необходимо дополнить задачу следующими граничными условиями.

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

Теперь мы можем перейти к программированию метода решения задачи (1-4) на основе явных схем первого порядка аппроксимации.

3. Программная реализация

Ниже приведена примерная реализация исключительно самого метода. Всё, что касается визуализации, решено отдельно.

Первый блок кода, это определение граничных условий.

#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);
}

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

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];
                }
        }
}

Здесь реализован вариант с устойчивой схемой -- разностным оператором D'x

4. Результаты

На следующей картинке нарисованы значения функции u(x,t) при a=0.5 и t=0, t=5.0 и t=10.0.

Видно, что волна не столько бежит, сколько размазывается в пространстве.

Теперь заменим основной цикл, на цикл, реализующий абсолютно неустойчивую схему, определенную оператором Dx.

    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);
    }

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

Этот результат совсем не похож на бегущую вправо ступеньку.

5. Заключение

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

Answer 2

Раз уж никто не желает разбираться с этим вопросом, то осмелюсь высказать свои соображения.

В вопросе не хватает части кода, а именно:

1. Не заданы значения всем элементам массивов wx, wt и wht. Когда объявляется массив, ему выделяется место в памяти и, если после объявления не задать значения элементам массива, то при обращении к ним получаешь то, что там оказалось. Отсюда и берутся эти "страшные" значения. ))
2. Где объявляется переменная T ??
3. Где объявляется переменная a ??

Несколько измененный код:

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]);
    }
}