1. Задание начального шаг - h, начальных значений xо,y10,…,yN0 и точности вычисления - ε.
2. В подпрограмме-процедуре задаём вид системы дифференциальных уравнений
3. В подпрограмме-функции задаём вид правой части уравнений
4. С помощью пяти циклов с управляющей переменной J=1,N вычисляем коэффициенты по формулам (3)-(7).
5. В последнем цикле находим решение системы дифференциальных уравнений по формуле (8) и погрешность по формуле (9).
6. Проверка выполнение условий (10) и (11). Если первое условие не выполняется то h:=h/2 и переходим к п.3
7. Если выполняются оба условия, то значение xi+1=xi+h и Yj(i+1) выводим на экран.
8. Если второе условие не выполняется, то h:=h+h и переходим к п.3
9. Вывести результаты вычислений в новом окне.
ГЛАВА II. ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРЕМЕНТ
Постановка задачи
Ставится задача составить программу решения системы дифференциальных уравнений на примере:
(12)
Начальные условия:
y(1)=2; y(2)=1. (13)
Требуется найти решение системы дифференциальных уравнений (12) с начальными условиями (13) методом Рунге-Кутта-Мерсона.
Анализ результатов
Система обыкновенных дифференциальных уравнений была решена методом Рунге-Кутта-Мерсона (см. Приложение 1).
Для проверки результата выполнения программы данная система была решена в MathCad(см. Приложение 2).
Для анализа результатов построим графики данных систем:
Сравнение результатовпоказывает, что они обеспечивают примерно одинаковое решение.
В большинстве случаев метод Рунге-Кутта-Мерсона даёт более точный результат (погрешность 
). Кроме того, хотя он громоздок в реализации, но быстрая сходимость метода компенсирует увеличение числа вспомогательных операций и, резко уменьшает вероятность числовой неустойчивости.
ЗАКЛЮЧЕНИЕ
В данной курсовой работе реализована средствами языка программирования Delphi программа, позволяющая решить систему обыкновенных дифференциальных уравнений методом Рунге-Кутта-Мерсона.
Также было проверенно решение данной системы в MathCad и проанализированы результаты.
Из анализа результатов вычисления можно сделать вывод о большей точности вычисления по методу Рунге-Кутта-Мерсона.
СПИСОК ЛИТЕРАТУРЫ
1. Н. Бахвалов, Н. Жидков, Г. Кобельков. Численные методы. М., 2002, 632 с.
2. Н. Калиткин. Численные методы. М., 1972,
3. А. Самарский. Введение в численные методы. М.,, 270с.
4. Ю. Тарасевич. Численные методы на Mathcad’е. Астрахань, 2000, 70с.
5. М. Лапчик, М. Рагулина, Е. Хеннер. Численные методы.М., 2004, 384с.
ПРИЛОЖЕНИЕ
Приложение 1
unit Unit1;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, ExtCtrls;
type
TForm1 = class(TForm)
Label1: TLabel;
Image1: TImage;
Label3: TLabel;
Edit1: TEdit;
Label4: TLabel;
Edit2: TEdit;
Label5: TLabel;
Label6: TLabel;
Edit3: TEdit;
Edit4: TEdit;
Button1: TButton;
Label2: TLabel;
Edit5: TEdit;
procedure Button1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
Const n=2;
Type
mas=Array[1..4] of String[60];
var
Form1: TForm1; i,s,p:integer;
h,v,x,e1,e2,e3: real;
y,w,k,f,e,a,c,d:array[1..10] of real;
j,d2:integer;
k2:byte;
implementation
uses Unit2;
{$R *.dfm}
function f1 (x1,y1,y2:real):real;
begin
f1:=y1+y2-x1*x1+x1-2;
end;
function f2(x1,y1,y2:real):real;
begin
f2:=-2*y1+4*y2+2*x1*x1-4*x1-7;
end;
procedure ur;
begin
f[1]:=f1(x,y[1],y[2]);
f[2]:=f2(x,y[1],y[2]);
end;
procedure TForm1.Button1Click(Sender: TObject);
begin
h:=strtofloat(edit1.Text);
e1:=strtofloat(edit2.Text);
x:=strtofloat(edit5.Text);
w[1]:=strtofloat(edit3.Text);
w[2]:=strtofloat(edit4.Text);
k2:=0;
e3:=0;
ur;
d2:=0;
for j:=1 to n do
begin
a[j]:=f[j]*H;
y[j]:=W[j]+a[j]/3;
end;
x:=x+h/3;
ur;
for j:=1 to n do
begin
y[j]:=W[j]+(a[j]+f[j]*H)/6;
end;
ur;
for j:=1 to n do
begin
c[j]:=f[j]*H;
y[j]:=W[j]+a[j]/8+0.375*c[j];
end;
x:=x+h/6;
ur;
for j:=1 to n do
begin
d[j]:=f[j]*H;
y[j]:=W[j]+a[j]/2-1.5*c[j]+2*d[j];
end;
x:=x+h/2;
ur;
for j:=1 to n do
begin
e[j]:=f[j]*H;
y[j]:=W[j]+(a[j]+4*d[j]+e[j])/6;
e2:=abs(-2*a[j]+9*c[j]-8*d[j]+e[j])/30;
if e2<=e1 then
if e2<e1/20 then d2:=d2+1 else
e3:=0;
end;
if e3<>0 then begin
x:=x-h;
for j:=1 to n do begin
y[j]:=W[j];
end;
H:=H/2;
end
else k2:=1;
if d2=n then H:=H+H;
form2.Show;
form2.edit1.text:=floattostr(y[1]);
form2.edit2.text:=floattostr(y[2]);
end;
end.
Приложение 2
|
|
|
|
|
|