Численные методы решения краевых задач для дифференциальных уравнений

13

Введение

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

Также данная курсовая работа включает в себя: описание метода, применение метода к конкретной задаче (анализ), код программы решения вышеперечисленных методов на языке программирования Borland C++ Builder 6.

Описание метода

Наша задача состоит в том, чтобы рассмотреть и реализовать в программном виде методы решения первой краевой задачи конечно-разностным методом и методом прогонки.

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

Теперь отдельно по каждому методу.

1. Конечно-разностный метод.

Сначала рассмотрим конечно-разностный метод решения краевой задачи.

Пусть дано уравнение вида

y''=b(x)y'(x)+c(x)y(x)+d(x) (1)

определенного не отрезке x [a,b], граничные условия y(a)= , y(a)= (2). Форма разделения из [a,b] используется в точках a=x0 < x1< … < xn =b, с шагом h=(b-a)/n и xk =a+kh, k=0,1,…,n.

Разностная схема данного дифференциального уравнения

(3)

(4)

Стартовые производные y(xk) мы заменяем на yk , в результате получаем отношение:

(5)

На следующем шаге мы избавляемся от O(h?) в (5) и вводим , и ; получаем дифференциальное уравнение:

, (6)

которое считает числовую аппроксимацию для дифференциального уравнения (1).

Преобразуем (6) следующим образом:

(7)

где k=1,2,…,n-1, y0 = , yn = .

Запишем (7) в матричном виде

где , .

2. Метод прогонки.

Общим классическим методом численного решения краевых задач с любым типом краевых условий является метод прогонки. Рассмотрим разностный подход к решению краевой задачи методом прогонки.

Также рассмотрим уравнение (1) на отрезке x[0,1] с граничным условием y(0)= , y(1)= (для 1-й краевой задачи).

Разностная схема данного дифференциального уравнения

(8)

преобразуем (8) следующим образом

(9)

то есть приводится (9) к виду

Ak yk-1 - Ck yk + Bk yk+1 = -fk (10)

(11)

(12)

Для 1-й краевой задачи ; .

Мы будем использовать решение из (10)-(12) в форме

(13)

где , .

Применение метода к конкретной задаче (анализ)

Дано дифференциальное уравнение вида 4xy'' + 2y' + y = 0 (14) , y(1)=1.3817 , x[1,2]. Точное решение этого уравнения имеет вид .

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

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

и решаем эту систему рассмотренным ранее методом Гаусса.

Теперь рассмотрим метод прогонки.

Уравнение (14) запишем в виде , преобразуем его к виду , где Ak = , Bk =, C==>

Ak yk-1 - Ck yk + Bk yk+1=0 и, используя уравнение (13) находим решение дифференциального уравнения.

Листинг программы:

//---------------------------------------------------------------------------

#include <vcl.h>

#pragma hdrstop

#include "UBound.h"

#include <math.h>

//---------------------------------------------------------------------------

#pragma package(smart_init)

#pragma resource "*.dfm"

TForm1 *Form1;

int n=10;

int x1=1,x2=2;

double h=(x2-x1)/(double)n;

//---------------------------------------------------------------------------

__fastcall TForm1::TForm1(TComponent* Owner)

: TForm(Owner)

{

}

//---------------------------------------------------------------------------

double ex_sol(double x)

{

return sin(sqrt(x))+cos(sqrt(x));

}

void __fastcall TForm1::Button_BashClick(TObject *Sender)

{

Label_TypeMethod->Caption="Конечно-разностный";

double x,s,p,q,par_1,par_n;

int i,j,i1,j1;

double A[100][100];

double y[100];

x=x1+h;

par_1=(h/2.0*(-1.0/(2.0*(x)))+1)*ex_sol(1);

par_n=(-h/2.0*(-1.0/(2.0*(2-2*h)))+1)*ex_sol(2);

for(i=0;i<n-1;i++)

{

q=(-1.0/(4.0*x));

p=(-1.0/(2.0*x));

for(j=0;j<n-1;j++)

{

if(i==j)

{

A[i][j]=(2+h*h*q);

}

else

{

if(j==(i+1))

{

A[i][j]=(h*p/2-1);

}

else

{

if(i==j+1)

A[i][j]=(-h*p/2-1);

else

A[i][j]=0;

}

}

}

if(i==0)

{

A[i][j]=par_1;

}

else

{

if(i==n-2)

A[i][j]=par_n;

else

A[i][j]=0;

}

x=x+h;

}

for(i=0;i<n-1;i++)

{

j=i;

s=A[i][j];

for(j=i;j<=n-1;j++)

{

A[i][j]=A[i][j]/s;

}

for(i1=i+1;i1<n-1;i1++)

{

s=A[i1][i];

for(j1=i;j1<n;j1++)

A[i1][j1]=A[i1][j1]-(A[i][j1]*s);

}

}

y[n-2]=A[n-2][n-1]/A[n-2][n-2];

for(i=n-3;i>=0;i--)

{

s=A[i][n-1];

for(j=i+1;j<=n-2;j++)

{

s=s-A[i][j]*y[j];

}

y[i]=s/A[i][i];

}

x=x1+h;

for (int m = 0; m < n-1; m++)

{

StringGrid1->Cells[0][m]=FloatToStr(x);

StringGrid1->Cells[1][m]=FloatToStr(y[m]);

StringGrid1->Cells[2][m]=FloatToStr(ex_sol(x));

StringGrid1->Cells[3][m]=FloatToStr(fabs(y[m]-ex_sol(x)));

x=x+h;

}

}

//---------------------------------------------------------------------------

void __fastcall TForm1::Button_ExitClick(TObject *Sender)

{

Close();

}

//---------------------------------------------------------------------------

void __fastcall TForm1::Button_MolClick(TObject *Sender)

{

Label_TypeMethod->Caption="Метод прогонки";

double x;

int i;

i = 0;

double alfa[100];

double beta[100];

double A,B,C,F;

alfa[0] = 0;

beta[0] = 1.38177;

x=x1;

for (i = 0; i < n - 1; i++)

{

x = x + h;

A=-h/2.0*(-1.0/(2.0*x))-1.0;

B=h/2.0*(-1.0/(2.0*x))-1.0;

C=-2.0-h*h*(-1.0/(4.0*x));

F=0;

alfa[i+1]=B/(C-alfa[i]*A);

beta[i+1]=(F+beta[i]*A)/(C-alfa[i]*A);

}

double y[100];

y[0] = ex_sol(x2);

x = x1+h;

for (i = 1; i < n; i++)

{

y[i] = alfa[n - i] * y[i - 1] + beta[n - i];

}

for (i = 0; i <=n; i++)

{

StringGrid1->Cells[0][i]=FloatToStr(x);

StringGrid1->Cells[1][i]=FloatToStr(y[n-i-1]);

StringGrid1->Cells[2][i]=FloatToStr(ex_sol(x));

StringGrid1->Cells[3][i]=FloatToStr(fabs(y[n-i-1]-ex_sol(x)));

x=x+h;

}

}

//---------------------------------------------------------------------------

void __fastcall TForm1::RadioButton1Click(TObject *Sender)

{

Button_Bash->Visible=true;

Button_Mol->Visible=false;

}

//---------------------------------------------------------------------------

void __fastcall TForm1::RadioButton2Click(TObject *Sender)

{

Button_Bash->Visible=false;

Button_Mol->Visible=true;

}

//---------------------------------------------------------------------------

void __fastcall TForm1::Button_ExitClick(TObject *Sender)

{

Close();

}

//---------------------------------------------------------------------------

void __fastcall TForm1::Button_MolClick(TObject *Sender)

{

Label_TypeMethod->Caption="Метод прогонки";

double x;

int i;

i = 0;

double alfa[100];

double beta[100];

double A,B,C,F;

alfa[0] = 0;

beta[0] = 1.38177;

x=x1;

for (i = 0; i < n - 1; i++)

{

x = x + h;

A=-h/2.0*(-1.0/(2.0*x))-1.0;

B=h/2.0*(-1.0/(2.0*x))-1.0;

C=-2.0-h*h*(-1.0/(4.0*x));

F=0;

alfa[i+1]=B/(C-alfa[i]*A);

beta[i+1]=(F+beta[i]*A)/(C-alfa[i]*A);

}

double y[100];

y[0] = ex_sol(x2);

x = x1+h;

for (i = 1; i < n; i++)

{

y[i] = alfa[n - i] * y[i - 1] + beta[n - i];

}

for (i = 0; i <=n; i++)

{

StringGrid1->Cells[0][i]=FloatToStr(x);

StringGrid1->Cells[1][i]=FloatToStr(y[n-i-1]);

StringGrid1->Cells[2][i]=FloatToStr(ex_sol(x));

StringGrid1->Cells[3][i]=FloatToStr(fabs(y[n-i-1]-ex_sol(x)));

x=x+h;

}

}

//---------------------------------------------------------------------------

void __fastcall TForm1::RadioButton1Click(TObject *Sender)

{

Button_Bash->Visible=true;

Button_Mol->Visible=false;

}

//---------------------------------------------------------------------------

void __fastcall TForm1::RadioButton2Click(TObject *Sender)

{

Button_Bash->Visible=false;

Button_Mol->Visible=true;

}

//---------------------------------------------------------------------------

Результаты расчета:

1. Конечно - разностный метод

x	exact solution	y	error
1,1	1,36543385512386	1,3655241231777	9,02680538475652E-5
1,2	1,3467827931773	1,34696264285001	0,0001798496727159
1,3	1,32614183910416	1,32640987061151	0,000268031507356567
1,4	1,30377563353274	1,32640987061151	0,000354381017934105
1,5	1,27990532272839	1,28034396043135	0,000438637702957411
1,6	1,25471820307239	1,25523885241359	0,000520649341198885
1,7	1,2283747392501	1,2289750725619	0,000600333311800857
1,8	1,20101378595756	1,20169143839785	0,000677652440289624
1,9	1,17275654886199	1,17350914840487	0,000752599542883688

2. Метод прогонки

x	exact solution	y	error
1,1	1,36543385512386	1,36542117841512	1,26767087399772E-5
1,2	1,3467827931773	1,34676470183753	1,80913397702582E-5
1,3	1,32614183910416	1,32612121470833	2,06243958234906E-5
1,4	1,30377563353274	1,303754611594	2,10219387379395E-5
1,5	1,27990532272839	1,27988551277913	1,98099492652827E-5
1,6	1,25471820307239	1,25470083439625	1,73686761431049E-5
1,7	1,2283747392501	1,22836076021839	1,39790317067597E-5
1,8	1,20101378595756	1,20100393346566	9,8524918956545E-6
1,9	1,17275654886199	1,17275139794799	5,15091399884513E-6