Барон Яков. Пример отчета по практическому заданию 1

1. Постановка задачи.Используя схему бегущего счета и итерационные методы, решить задачу:2. Метод решения.Введем одномерные пространственную и временную сетки в областигде:Перепишем исходное уравнение в виде:На введенной разностной сетке будем рассматривать сеточную функциюдифференциальные операторы на конечно-разностные:ЗаменимДля граничных условий будем иметь:В данной задаче будем использовать трёхточечный шаблон с порядком аппроксимации:Пусть:m+1Решение данного уравнения будем искать методомпоследовательных приближений. В качестве первогоприближения возьмем значение функции в узле.Для n+1 приближения будем иметь:Вычисления прекращаются при условии, чтозаданное число (точность приближенного решения).mnn+1., гденаперед13. Характеристики уравнения.Пусть.

Характеристики данного квазилинейного уравнения переносаопределяются соотношениями:Для того, чтобы решение не претерпевало разрыв, характеристики не должныпересекаться. Учитывая начальные условия, получим::Графики характеристик представлены на рисунке:Характеристики пересекаются примерно при.4. Устойчивость разностной схемы.Выберем произвольную внутреннюю точку области, где рассматривается поставленнаязадача, и применим метод «замороженных» коэффициентов в данной точке.2деИщем решение уравнения в виде (спектральное разложение):подстановки в вышеуказанное выражение, получим выражение для :. ПослеИз данного выражения нетрудно получить оценку:Таким образом, выполнено необходимое условие Неймана.Выполним проверку достаточного условия.

Выпишем наше разностное уравнение внесколько иной форме:В силу последней оценки оказывается выполненным достаточное условие устойчивости(критерий Куранта).35. Численное решение.Ниже представлен код программы, реализующей численное решение исходной задачи.Программа выполнена на языке C. Графическая визуализация получена с помощьюпрограммного пакета ROOT.#include "stdio.h"#include <stdlib.h>#include "iostream"#include "clocale"#include "math.h"#include "iomanip"#include "fstream"#define pi 3.1415926535using namespace std;int Nx, Nt;double delta = 0.0001;double **U;double dx, dt;void base() {U=(double**)malloc(Nt*sizeof(double*));for (int i=0; i<=Nt; i++) {U[i]=(double*)malloc(Nx*sizeof(double));}}void bounds() {for(int n=0; n<=Nx; n++) {U[0][n]=sin(pi*(n*dx));}for(int m=0; m<=Nt; m++) {U[m][0]=exp(-m*dt)-1;}}double f(double u) {return (atan(exp(u)));}double f(int i, int j) {return f(U[i][j]);}double df(double u) {return exp(u)/(1+exp(2*u));}double equation(double u, int i, int j) {return ((u-U[i][j+1])/dt+(f(u)-f(U[i+1][j]))/dx);}double dif_eq (double u) {return 1/dt + df(u)/dx;}double solution(int i, int j) {double u_prev=1, u=0;while (abs(u-u_prev)>delta) {u_prev = u;u = u_prev - equation(u_prev, i ,j)/dif_eq (u_prev);}4return u;}void maker() {for(int i = 0; i < Nt-1; i++) {for(int j = 0; j < Nx-1; j++) {U[i+1][j+1]=solution(i, j);}}}int OMM1 (void) {Nx = 100; Nt = 100;dx = (1./Nx), dt = (1.8/Nt);base();bounds();maker();return 0;}56.

Результат численного решения.void picture() {Double_t x=0, t=0;TCanvas *c = new TCanvas("c","u", 0, 0, 600, 400);TGraph2D *a = new TGraph2D();int N=1;for (int i = 0; i <= Nt; i++) {x=0;for (int j = 0; j <= Nx; j++) {x+=dx;a->SetPoint(N, t, x, U[i][j]);N++;}t+=dt;}gStyle->SetPalette(1);a->GetXaxis()->SetTitle("x");a->GetYaxis()->SetTitle("t");a->GetZaxis()->SetTitle("U");a->Draw("surf1");gStyle->SetPalette(1);TCanvas *c1 = new TCanvas("c1","u1", 0, 0, 600, 400);TGraph *a1 = new TGraph();N=1;for (int j = 0; j <= Nx; j++) {x=j*dx;a1->SetPoint(N,x,U[0][j]);N++;}a1->GetXaxis()->SetTitle("x");a1->GetYaxis()->SetTitle("U");a1->Draw();}6.