Многошаговые методы решения дифференциальных уравнений

12

Введение

Данная курсовая работа включает в себя три многошаговых метода решения обыкновенного дифференциального уравнения.

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

1. Явный метод Адамса - Башфорта.

2. Неявный метод Адамса - Моултона.

3. Метод прогноза и коррекции.

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

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

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

Рассмотрим задачу Коши:

y' = f(x,y) , a? x?b (1)

y(a) =y₀ (2)

Точность вычислений в отличие от одношаговых методов можно увеличить, если использовать информацию о нескольких предыдущих точках x_k , x_k-1 , …

Большой и важный класс многошаговых методов возникает на основе следующего подхода. Если подставить в формулу (1) точное решение y(x) и проинтегрировать это уравнение на отрезке (x_k , x_k+1 ), то получим:

, (3)

где в последнем члене предполагаем, что p(x) - полином, аппроксимирующий f(x,y(x)) . Чтобы построить этот полином, предположим, что y_k , y_k-1 , … , y_k-n - приближения к решению в точках x_k , x_k-1 , … , x_k-n . Мы по-прежнему считаем, что узлы расположены равномерно с шагом h. Тогда f_i = f(x_i , y_i ) (i = k, k-1, …, k-n) есть приближения к f(x,y(x)) в точках x_k , x_k-1 , … , x_k-n , и мы в качестве p возьмём полином для набора данных (x_i , f_i ) (i = k, k-1, …, k-n) . Таким образом, p - полином степени n , удовлетворяющий условиям p(x_i ) = f_i , (i = k, k-1, …, k-n). В принципе, можно проинтегрировать этот полином явно, что ведёт к следующему методу:

 (4)

В простейшем случае, когда n=0 , полином p есть константа, равная f_k и (4) превращается в обычный метод Эйлера. Если n=1 , то p есть линейная функция, проходящая через точки (x_{k -1} , f_{k -1} ) и (x_k , f_k ) ,т.е.

.

Интегрируя этот полином от x_k до x_{k +1,} получаем следующий метод:

, (5)

который является двухшаговым, поскольку использует информацию в двух точках x_k и x_{k -1} .

Аналогично, если n=2,то p есть квадратичный полином, интерполирующий данные (x_{k -2} , f_{k -2} ) , (x_{k -1} , f_{k -1} ) и (x_k , f_k ) , а метод имеет вид:

(6)

Если n=3 , то интерполяционный полином является кубическим, а соответствующий метод определяется формулой:

(7)

Методы (6) и (7), соответственно называются трёхшаговым и четырёхшаговым.

Формулы (5)-(7) известны как явные методы Адамса (Адамса-Башфорта), т.к. они для нахождения y_k+1 не требуют решения никаких уравнений. Метод (5) имеет второй порядок точности, поэтому его называют методом второго порядка. Аналогично, методы (6) и (7) называют соответственно методами Адамса-Башфорта третьего и четвёртого порядков.

Методы Адамса-Башфорта используют уже сосчитанные значения в точке x_k и в предыдущих точках. В принципе, при построении интерполяционного полинома мы можем использовать и точки x_k+1 , x_k+2 и т.д. Простейший случай при этом состоит в использовании точек x_k+1 , x_k , … , x_k-n и построении интерполяционного полинома степени n+1 , удовлетворяющего условиям p(x_i )= f_i (i = k+1, k, …, k-n) . При этом возникает класс методов, известных как неявные методы Адамса (Адамса-Моултона).

Если n=0, то p - линейная функция, проходящая через точки (x_k , f_k ) и (x_k+1 , f_k+1 ) , и соответствующий метод

(8)

является методом Адамса-Моултона второго порядка.

Если n=2, то p - кубический полином, построенный по точкам (x_k+1 , f_k+1 ) , (x_k , f_k ), (x_{k -1} , f_{k -1} ) и (x_{k -2} , f_{k -2} ), и соответствующий метод

(9)

является методом Адамса-Моултона четвёртого порядка.

Заметим теперь, что в формулах (8) и (9) значение f_k+1 неизвестно. Дело в том, что для вычисления f(x_k+1 , y_k+1 )= f_k+1 нужно знать значение y_k+1 , которое само пока является неизвестным. Следовательно, методы Адамса-Моултона определяют y_k+1 только неявно. Так, например, соотношение (8) действительно является уравнением относительно неизвестного значения y_k+1. То же самое справедливо и относительно (9). В силу этого методы Адамса-Моултона называются неявными. В то же время методы Адамса-Башфорта называют явными, поскольку они для нахождения значения y_k+1 не требуют решения никаких уравнений.

Теперь рассмотрим метод прогноза и коррекции. Он также относится задачам класса Коши, а именно к численным решениям многошаговыми методами.

Рассмотрим 2-х шаговый метод прогноза и коррекции. Пусть дано ДУ для которого известно значение функции в двух соседних узлах сетки

, (xk-1, yk-1) и (xk, yk).

Сначала строится прогноз значения по формуле Адамса-Башфорта

.

Затем это значение корректируется по более точной формуле.

Для того чтобы начать расчет методом прогноза и коррекции, необходимо знать значения функции в двух первых узлах сетки - x₀ и x₁ - (x0, y0) и (x1, y1). Обычно значение в узле x₁ определяется каким-либо одношаговым методом.

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

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

Дано дифференциальное уравнение вида

4xy'' + 2y' + y = 0, x[1,2], где y(1)=1,38177, y'(1)=-0,15058.

Точное решение этого уравнения имеет вид .

Задача состоит в том, чтобы решить это дифференциальное уравнение, применив многошаговые методы решения ОДУ: явный метод Адамса-Башфорта, неявный метод Адамса-Моултона, метод прогноза и коррекции.

Как было описано выше данное уравнение можно решить как методом второго порядка, так и методами третьего и четвертого порядков. Но я применил методы второго порядка во всех трех случаях (в случае метода Адамса-Башфорта применил формулу (5), в Адамсе-Моултоне применил формулу (8), а в методе прогноза и коррекции прогноз строится по формуле (5), а коррекция по формуле (8)).

Изучив методы и применив их к данному дифференциальному уравнению, приходим к такому выводу: при решении данного уравнения тремя известными способами при шаге h=0,1 все три метода дают хорошую точность.

Достоинством метода решения ОДУ заключается в том, что в каждой точке рассчитывается только одно значение функции f(x,y) и они являются весьма экономичными по затратам времени вычислений при достаточной точности.

К недостаткам можно отнести невозможность старта многошагового метода из единственной начальной точки, так как для вычислений по k-шаговой формуле необходимо значение значения функции в k узлах. Поэтому приходится (k-1) решение в первых узлах x₁, x₂, …, x_k-1 получать с помощью какого-либо одношагового метода (я применил метод Рунге-Кутте 2-го порядка). Другим недостатком является невозможность изменения шага в процессе решения, что легко реализуется в одношаговых методах.

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

#include <vcl.h>

#pragma hdrstop

#include "UMoreStep.h"

#include <math.h>

#pragma package(smart_init)

#pragma resource "*.dfm"

TForm1 *Form1;

int n=10;

double x1=1,x2=2;

double x[100],y[100],p_y[100];

double h=(x2-x1)/n;

__fastcall TForm1::TForm1(TComponent* Owner)

: TForm(Owner)

{

}

double dif_ur(double x,double y,double p_y)

{

return ((-2*p_y-y)/(4*x));

}

void __fastcall TForm1::Button_BashClick(TObject *Sender)

{

Label_TypeMethod->Caption="Aaain-Aaooi?o";

y[0] = 1.38177

x[0] = x1;

p_y[0] = -0.15058;

x[1] = x[0] + h;

p_y[1]=p_y[0]+(h/2)*(dif_ur(x[0],y[0],p_y[0])+dif_ur(x[0]+h,y[0],p_y[0]+h*dif_ur(x[0],y[0],p_y[0])));

y[1]=y[0]+(h/2)*(p_y[0]+p_y[1]);

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

{

x[i+1]=x[i]+h;

p_y[i+1]=p_y[i]+h/2*(3*dif_ur(x[i],y[i],p_y[i])-dif_ur(x[i-1],y[i-1],p_y[i-1]));

y[i+1]=y[i]+h/2*(3*p_y[i]-p_y[i-1]);

}

for (int j = 0; j <= n; j++)

{

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

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

StringGrid1->Cells[2][j]=FloatToStr(sin(sqrt(x[j]))+cos(sqrt(x[j])));

StringGrid1->Cells[3][j]=FloatToStr(fabs(StrToFloat(StringGrid1->Cells[2][j])-y[j]));

}

}

void __fastcall TForm1::Button_ExitClick(TObject *Sender)

{

Close();

}

void __fastcall TForm1::Button_MolClick(TObject *Sender)

{

Label_TypeMethod->Caption="Aaain-Iioeoii";

y[0] = 1.38177;

x[0] = x1;

p_y[0] = -0.15058;

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

{

x[i+1]=x[i]+h;

p_y[i+1] = p_y[i] + (h * dif_ur(x[i], y[i], p_y[i]) + h * dif_ur(x[i] + h, y[i] + h * p_y[i], p_y[i] + h * dif_ur(x[i], y[i], p_y[i]))) / 2;

y[i+1] = y[i] + h * (p_y[i] + h * dif_ur(x[i], y[i], p_y[i]) / 2);

p_y[i+1] = (p_y[i] + h * (dif_ur(x[i] + h, y[i+1], p_y[i+1]) + dif_ur(x[i], y[i], p_y[i])) / 2);

y[i+1] = y[i] + h * (p_y[i+1] + p_y[i]) / 2;

}

for (int j = 0; j <= n; j++)

{

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

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

StringGrid1->Cells[2][j]=FloatToStr(sin(sqrt(x[j]))+cos(sqrt(x[j])));

StringGrid1->Cells[3][j]=FloatToStr(fabs(StrToFloat(StringGrid1->Cells[2][j])-y[j]));

}

}

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

void __fastcall TForm1::RadioButton1Click(TObject *Sender)

{

Button_Bash->Visible=true;

Button_Mol->Visible=false;

Button_Cor->Visible=false;

}

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

void __fastcall TForm1::RadioButton2Click(TObject *Sender)

{

Button_Bash->Visible=false;

Button_Mol->Visible=true;

Button_Cor->Visible=false;

}

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

void __fastcall TForm1::RadioButton3Click(TObject *Sender)

{

Button_Bash->Visible=false;

Button_Mol->Visible=false;

Button_Cor->Visible=true;

}

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

void __fastcall TForm1::Button_CorClick(TObject *Sender)

{

Label_TypeMethod->Caption="I?aaeeoi?-ei??aeoi?";

y[0] = 1.38177;

x[0] = x1;

p_y[0] = -0.15058;

x[1] = x[0] + h;

p_y[1]=p_y[0]+(h/2)*(dif_ur(x[0],y[0],p_y[0])+dif_ur(x[0]+h,y[0],p_y[0]+h*dif_ur(x[0],y[0],p_y[0])));

y[1]=y[0]+(h/2)*(p_y[0]+p_y[1]);

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

x[i+1]=x[i]+h;

p_y[i+1] = p_y[i] + h * (3 * dif_ur(x[i],y[i],p_y[i]) - dif_ur(x[i-1],y[i-1],p_y[i-1])) / 2;

y[i+1] = y[i] + h * (3* p_y[i] - p_y[i-1]) / 2;

p_y[i+1] = p_y[i] + h * (dif_ur(x[i]+h, y[i+1], p_y[i+1]) + dif_ur(x[i], y[i], p_y[i]))/2;

y[i+1] = y[i] + h * (p_y[i+1] + p_y[i])/2;

}

for (int j = 0; j <= n; j++)

{

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

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

StringGrid1->Cells[2][j]=FloatToStr(sin(sqrt(x[j]))+cos(sqrt(x[j])));

StringGrid1->Cells[3][j]=FloatToStr(fabs(StrToFloat(StringGrid1->Cells[2][j])-y[j]));

}

}

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

Метод Адамса-Башфорта

x	exact solution	y	error
1	1,38177329067604	1,38177	3,29067604005725E-6
1,1	1,36543385512386	1,36545333494318	1,94798193216625E-5
1,2	1,3467827931773	1,34661933977273	0,000163453404572925
1,3	1,32614183910416	1,32588519852612	0,000256640578043141
1,4	1,30377563353274	1,30346514735694	0,000310486175804584
1,5	1,27990532272839	1,27957128853577	0,000334034192624262
1,6	1,25471820307239	1,25438392699673	0,000334276075660152
1,7	1,2283747392501	1,22805836929195	0,000316369958149775
1,8	1,20101378595756	1,20072959517654	0,000284190781024675
1,9	1,17275654886199	1,17251584901979	0,000240699842204883
2	1,14370964075811	1,1435214460441	0,000188194714010557

Метод Адамса-Моултона

x	exact solution	y	error
1	1,38177329067604	1,38177	3,29067604005725E-6
1,1	1,36543385512386	1,36546037053267	2,65154088104422E-5
1,2	1,3467827931773	1,34682666087201	4,38676947109723E-5
1,3	1,32614183910416	1,32619402286957	5,21837654078901E-5
1,4	1,30377563353274	1,30382943566818	5,38021354416577E-5
1,5	1,27990532272839	1,27995569425606	5,03715276659314E-5
1,6	1,25471820307239	1,25476128734688	4,30842744907492E-5
1,7	1,2283747392501	1,22840756096378	3,28217136845617E-5
1,8	1,20101378595756	1,20103403374354	2,02477859798972E-5
1,9	1,17275654886199	1,17276241989401	5,87103201764609E-6
2	1,14370964075811	1,14369972745694	9,91330116907119E-6

Метод прогноза и коррекции

x	exact solution	y	error
1	1,38177329067604	1,38177	3,29067604005725E-6
1,1	1,36543385512386	1,36545333494318	1,94798193216625E-5
1,2	1,3467827931773	1,34680554703575	2,27538584458425E-5
1,3	1,32614183910416	1,32615890975456	1,70706503960812E-5
1,4	1,30377563353274	1,30378054845518	4,91492243921124E-6
1,5	1,27990532272839	1,27989334139954	1,19813288543E-5
1,6	1,25471820307239	1,25468582407744	3,23789949485528E-5
1,7	1,2283747392501	1,22831936753384	5,53717162598104E-5
1,8	1,20101378595756	1,20093350217706	8,02837804958173E-5
1,9	1,17275654886199	1,17264994579682	0,000106603065173971
2	1,14370964075811	1,1435757049175	0,000133935840614283