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

Обновил содержание вопроса,так как некоторые моменты я сам понял и исправил. И так есть следующие переменные:
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Если я напишу вместо 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;xSeries["Функция Хэвисайда"]->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 for(int j = 0;j ///-----------------Вычисление основых формул с разностной схемой
for(int j = 0;j //----Для записи в простой текстовый файл String^ fileName = "results.txt"; StreamWriter^ sw = gcnew StreamWriter(fileName); for(int j = 0;jWrite("{0} ",wht[j][i]); } sw->WriteLine(); } sw->Close();
//Тестовые функции для отображения drawfunc(T);// Вызов функции для рисования
///---------------Построение графика
for(int j = 0;jSeries["Series2"]->BorderWidth=3; chart2->Series["Series2"]->Points->AddXY(i,wht[j][i]); } } } else MessageBox::Show("Ошибка!Не все поля заполнены!"); }

Ответ

Для того, чтобы ответить, почему не бежит волна на графике, нужно сначала разобраться с численным методом, а лишь потом с кодом на 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 #include #include
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 for(int j=0; j < N-1; j++){ u[j+1][0] = u1(xmin,j*ht); for(int i=1; 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Результат будет выглядеть примерно так

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