Программа для решения дифференциальных уравнений первого порядка методом Рунге-Кутта (стр. 4 из 6)
yi+1 = yi + (xi+1 - xi)f(xi, yi) = yi + hf(xi, yi).
Из полученного выражения видно, что вычисление интеграла по методу прямоугольников приводит к формуле Эйлера.
Вычислим интеграл по формуле трапеций:
yi+1 = yi +0,5h(f(xi, yi)+ f(xi+1, yi+1))
Из выражения видно, что оно совпадает с расчетной формулой усовершенствованного метода Эйлера-Коши. Для получения более точного решения дифференциального уравнения следует воспользоваться более точными методами вычисления интеграла. В методе Рунге-Кутта искомый интеграл представляется в виде следующей конечной суммы:
где Pi— некоторые числа, зависящие от q; Ki(h) — функции, зависящие от вида подынтегральной функции f(x,y) и шага интегрирования h, вычисляемые по следующим формулам:
K1(h) = hf(x, y);
K2(h) = hf(x + a2h, y + β21K1(h));
K3(h) = hf(x + a3h, y + β 31K1(h) + β 32K2(h));
Kn(h) = hf(x + aqh, , y + β q1K1(h) + ... + βq,q-1Kq-1(h)).
Значения p, α, β получают из соображений высокой точности вычислений. Формулы Рунге-Кутта третьего порядка (q= 3) имеют следующий вид:
K1=hf(xi, yi);
K2=hf(xi + 0,5h, yi+0,5 K1);
K3=hf(xi+h, yi+K1+2K2).
Наиболее часто используется метод Рунге-Кутта четвертого порядка, для которого расчетные формулы имеют следующий вид:
K1=hf(xi, yi);
K2=hf(xi + 0,5h, yi+0,5 K1);
K3=hf(xi+0,5h, yi+0,5K2).
K3=hf(xi+h, yi+K3).
Формулы Рунге-Кутта имеют погрешности порядка hq+1 . Погрешность метода Рунге-Кутта четвертого порядка имеет порядок h5
4.2 Описание программы ” РЕШЕНИЕ ОДУ “
Программа ”Решение ОДУ“ достаточно проста в использовании.
При запуске программы открывается главное окно программы (рис. 4), с установленными по умолчанию начальными условиями в полях ввода.
Назначение элементов ввода данных.
1. Поля X1, X2, Y(x1), Hпредназначены для ввода начального и конечного значений отрезка, на котором ищется решение дифференциального уравнения, значения функции при аргументе равном Х1 и величины шага дифференцирования;
2. В поле dY выводится формула дифференциального уравнения 1-й степени, выбранная для решения;
3. В поле dY(x1,y1) выводится значение производной в исходной точке.
Рис.4
Назначение элементов управления и контроля.
1. При нажатии кнопки EXAMPLE активируются “радиокнопки” выбора уравнений;
2. Щелчком “радиокнопки” выбирается соответствующее ей уравнение, вид формулы контролируется по её отображению в поле dY ;
3. Щелчком по кнопке ВЫЧИСЛИТЬ находятся приближенные решения выбранного дифференциального уравнения на заданном интервале;
4. Решения дифференциального уравнения в виде пар значений X - Y выводятся в поля X и Y; (рис. 5.)
По окончании вычислений активируются кнопка ГРАФИК и пункт меню ГРАФИК главного окна системы.
Рис.5
4.3 Назначение элементов графического окна программы
Вход в графическое окно осуществляется с помощью кнопок ГРАФИК наглавной формеили пункт меню ГРАФИК (рис. 6).
С помощью кнопки ВЫЧЕРТИТЬ на координатную плоскость выводится график функции – решение дифференциального уравнения на заданном интервале.
рис.6
4.4 Реакция программы при возникновении ошибок
При вводе пользователем ошибочных данных (отсутствии начальных условий, некорректных значений переменных) программа выдает сообщение об ошибке (рис.7 а, б) рис.7а. рис.7б.
рис 7а Рис.7б
Версия DELPHI
4.5 Перечень компонент DELPHI использованных впрограмме
В Form1 использованы компоненты:
- Edit1.text, Edit2.text, Edit3.text, Edit4.text – для ввода начальных условий дифференциального
уравнения
- Memo4.TMemo – для вывода формулы уравнения;
- Memo1.TMemo, Memo2.TMemo - для вывода результатов вычислений;
- Memo3.TMemo – для вывода значения производной в точке (Х0,Y0)
- ScrollBars ssVertical всвойствахMemo1.TMemo, Memo2.TMemo;
- Button1 “Вычислить”, Button2 “Очистить”, Button3 “График”, Button4 “Выход”,
Button5 “Example”, Button6 “UnExample”;
- Label1.TLabel - Label9.TLabel – дляотображенияназначениякомпонентов Memo иEdit;
- RadioGroup – для выбора вида уравнения;
- MainMenu ;
В Form2 использованы компоненты:
- MainMenu- для построения графика;
В Form3 использованы компоненты:
- Panel1.TPanel – для размещения информации о программе;
- Label1.TLabel – Label14.TLabel – для отображения информации о программе;
- Button1.TButton “OK” – для выхода из окна
5. ТЕХНИЧЕСКИЕ ХАРАКТЕРИСТИКИ И ТРЕБОВАНИЯ К ПО
Технические характеристики
Программа работает в среде операционных систем Windows 9х, NT.
Требования к ПО
Минимальные системные требования
aпроцессор Intel 486 с рабочей частотой 66 MHz и выше;
b) операционная система Windows 95, 98, NT 4.0, 2000, XP;
с) 16 Мбайт оперативной памяти (или более);
d) 3 Мбайт свободного пространства на жёстком диске.
6. ТЕКСТ ПРОГРАММЫ
Код программы
unitRKt;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls, CheckLst, ComCtrls, ExtCtrls,math, Menus;
type
TRKutta = class(TForm)
Label2: TLabel;
Label3: TLabel;
Label5: TLabel;
Label6: TLabel;
Edit1: TEdit;
Edit2: TEdit;
Edit3: TEdit;
Edit4: TEdit;
Memo1: TMemo;
Memo2: TMemo;
Button1: TButton;
Button2: TButton;
Label4: TLabel;
Label7: TLabel;
Button3: TButton;
Button4: TButton;
Label9: TLabel;
RadioGroup1: TRadioGroup;
Button5: TButton;
Memo3: TMemo;
Button6: TButton;
MainMenu1: TMainMenu;
N1: TMenuItem;
N3: TMenuItem;
N4: TMenuItem;
Example1: TMenuItem;
UnExample1: TMenuItem;
N5: TMenuItem;
N7: TMenuItem;
N2: TMenuItem;
N9: TMenuItem;
N8: TMenuItem;
Label1: TLabel;
Memo4: TMemo;
Label8: TLabel;
Label10: TLabel;
Label11: TLabel;
Label12: TLabel;
procedure Button1Click(Sender: TObject);
procedure Button2Click(Sender: TObject);
procedure Button3Click(Sender: TObject);
procedure Button4Click(Sender: TObject);
procedure RadioGroup1Click(Sender: TObject);
procedure Button5Click(Sender: TObject);
procedure Button6Click(Sender: TObject);
procedure N7Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
var
RKutta: TRKutta;
x1,x2,yc,xc,y,h:extended;
line_arr:integer; // размерность массивa
f:real; // значение функции при начальных условиях
zx:array of real;
zy:array of real;
implementation
uses RungeKutta, Spravka;
{$R *.DFM}
procedure TRKutta.Button1Click(Sender: TObject);
var k1,k2,k3,k4:extended;
t:integer;
begin
Memo1.Clear;
Memo2.Clear;
memo1.Enabled:=true;
memo2.Enabled:=true;
Memo3.lines.Add(''+floattostr(f));
Button1.Enabled:=False;
//Проверка возможности ввода начальных условий и инициализация переменных
try
// инструкции, которые могут вызвать исключение (ошибку)
x1:=StrToFloat(Edit1.Text); //инициализация(ввод) x1(начало отрезка)
xc:=StrToFloat(Edit1.Text); //передача начальных условий (х0)в графический модуль
x2:=StrToFloat(Edit2.Text); //инициализация(ввод) x2(конец отрезка)
y:=StrToFloat(Edit3.Text); //инициализация(ввод) Y(x1)
yc:=StrToFloat(Edit3.Text); //передача начальных условий (Y0)в графический модуль
h:=StrToFloat(Edit4.Text); //инициализация(ввод) H -величины шага вычислений
//определение размера массивов содержащего значения аргумента и функции
line_arr:=round(abs((x2-x1)/h))+1;
//Установка размера динамических массивов zx, zy
SetLength(zx, line_arr);
SetLength(zy, line_arr);
t:=0; // счётчик
while x1<x2
do
begin
zx[t]:=x1;
zy[t]:=y;
k1:=h*f;
if (y+(k1/2))=0 then begin showmessage ('Деление на 0!'+#13+' Измените h'); break;
end;
k2:=h*((y+(k1/2))-(2*(x1+(h/2)))/(y+(k1/2)));
if (y+(k2/2))=0 then begin showmessage ('Деление на 0!'+#13+' Измените h'); break;
end;
k3:=h*((y+(k2/2))-(2*(x1+(h/2)))/(y+(k2/2)));
if (y+k3)=0 then begin showmessage ('Деление на 0!'+#13+' Измените h');break;
end;
k4:=h*(y+k3-2*(x1+(h/2)))/(y+k3);
x1:=x1+h;
y:=y+(1/6)*(k1+2*k2+2*k3+k4);
t:=t+1;
Memo1.Lines.Add('x['+floattostr(t)+']='+floattostr(x1));
Memo2.Lines.Add('y['+floattostr(t)+']='+floattostr(y));
end;
except
onEConvertErrordo// невозможно преобразовать строку символов в число
begin
MessageDlg('Некорректные значения переменных',mtError,[mbOk],0);
exit;
end;
end;
{ | Подключение кнопок: |
| - Button3 - 'ГРАФИК': входа в графический модуль |
| - нопки N3 - 'EXAMPLE': образцовых функций | }
if ((y+(k1/2))=0) or ((y+(k2/2))=0) or((y+k3)=0) then
begin
Button3.Enabled:=False;
N3.Enabled:=False;
end
else
begin
Button3.Enabled:=True;
N3.Enabled:=True;
end;end;
//------------------------------------------------------------------------------
procedure TRKutta.Button2Click(Sender: TObject);
begin
memo1.Clear;
MEMO2.Clear;
MEMO3.Clear;
memo4.Clear;
edit1.Clear;
edit2.Clear;
edit3.Clear;
edit4.Clear;
end;
//------------------------------------------------------------------------------
{Процедура вывода окна ГРАФИК}
procedure TRKutta.Button3Click(Sender: TObject);
begin
Button3.enabled:=false;
N3.Enabled:=false;
Form2.ShowModal;
end;
//------------------------------------------------------------------------------
{Процедура выхода из программы}
procedure TRKutta.Button4Click(Sender: TObject);
begin
Close;
end;
//------------------------------------------------------------------------------
{Процедура выбора образцовой функции}
procedure TRKutta.RadioGroup1Click(Sender: TObject);
var x_rg,y_rg:extended;
begin
try
y_rg:=strtofloat(edit3.Text); //ввод Y(x1)
x_rg:=strtofloat(edit1.Text); //ввод X1
{------------------------------------------------------------------------------}
if RadioGroup1.ItemIndex=0 then
begin
if x_rg=0 then
begin
ShowMessage('введите X1 неравное 0');
RadioGroup1.Enabled:=False;
RadioGroup1.ItemIndex:=31;
end
else
begin
Memo3.clear;
Memo4.Clear;
Memo4.Lines.Add(' -(y+1)/x '); вывод формулы функции в окно "dY"
f:=-(y_rg+1)/x_rg; //вычисление значения dY