Реализация. 
Реализация метода Рунге-Кутты-Фельберга на языке программирования Delphi

Программа, реализующая адаптивную процедуру Рунге-Кутты 4 порядка. Написана программа на языке Delphi в среде CodeGear 2009.

//Класс для метода Рунге-Кутты-Фельдера.

TRunge = class.

private.

Erange: integer; //кол-вофункцийвсистеме.

Estart: TRealArr; //Начальные значения.

EFirst, ELast: double; //Начальное и конечное значение по оси t.

Eeps: double; //Точность.

Efun: TArrfun; //Массив функций.

EArr: TResArr; //Массив значений функции в точках, больше 0.

EArrNeg: TResArr; //Массив значений функции в точках, меньших 0.

EGetInd: integer; //Внутренняя переменная — для быстрого поиска нужного значения.

procedureSetRange (constValue: integer); virtual;//Установить количество функций в системе.

functionGetStart (Index: Integer): double; //Получить начальное значений для функции index.

procedureSetStart (Index: Integer; constValue: double);//Установить начальное значение value для функции index.

functionGetRes (X: double; index: integer): double; //Узнать значение функции index в точке x.

functionGetResInd (X: double): TRealArr; //Узнать значение всех функций в точке x.

functionGetCount: integer; //Узнать общее количество точек.

functionGetLeftInd: integer; //Узнать номер самой левой точки в массиве.

functionGetRightInd: integer; //Узнать номер самой правой точки.

functionValue (constx: real; consty: TRealArr; consth: real):TRealArr; // Вычисление значения в следующей точке методом РК-4.

functionValue5(constx: real; consty: TRealArr; consth: real):TRealArr;//Вычисление значения в следующей точке методом РК-5.

procedureSetEps (constValue: double); virtual; //Установить значение точности.

functioncalc (vlast: double; veps: double): TRes; overload; virtual; //Процедура вычисления значения в точке vlast с точностью veps.

public.

propertyLeftInd: integer read GetLeftInd; //Левыйиндекс.

propertyRightInd: integer read GetRightInd; // Правый индекс в массиве точек.

functionGetInd (v: integer; n: integer): double; //Узнать значение функции v в точке n.

propertyCount: integerreadGetCount; //Доступ к количеству точек в результате.

propertyArr[X: double]: TRealArrread GetResInd; default; //Доступ к значениям функции в точке x.

propertyArrInd[X: double; index: integer]: double read GetRes; //Доступ к значению конкретной функции index в точке x.

propertyrange: integer read Erangewrite SetRange; //Доступкколичествуфункций.

propertyStart[Index: integer]: double readGetStartwriteSetStart; //Доступ к начальным значениям.

property First: double readEfirstwriteEFirst; //Начальное значение по оси t.

property Last: double readELastwriteELast; //Конечное значение по оси t.

propertyeps: double readEEpswriteSetEps; //Точность вычисления.

property Fun: TarrFunreadEFunwriteEFun;// Доступ к массиву функций.

constructorCreate;virtual;

functionCalcat (vfirst, vlast, veps: double): boolean; //Вычислить значение на отрезке с точностью.

procedureProcessMessages;

functionCheckDelat (n: integer; f: TFunReal): double; //Проверить отклонение полученной функции от реальной.

end;

Класс TGraphicsRunge = class (TRunge) нужен для изображения графика.

//Класс для исходной задачи.

TProg1 = class (TRunge).

public.

constructorCreate; override; //В конструкторе описываем те входные данные, которые требуется использовать в конкретной задаче.

end;

//Два класса для тестовых примеров.

TProg2 = class (TRunge).

public.

constructor Create; override;

end;

TProg3 = class (TRunge).

public.

constructor Create; override;

end;

Конкретные задачи. В конструкторе определяются начальные условия.

Реализация метода Рунге-Кутты-Фельберга.

// Вычисление одного шага методом Р-К 4.

functionTRunge.Value (const x: real; const y: TRealArr; const h: real): TRealArr;

var i: integer; //счетчики.

k: arrayofTRealArr;

begin.

setlength (k, 4, Erange);

setLength (Result, Erange);

for i:=0 toErange — 1 do k[0, i]: =Efun[i](x, y);

for i:=0 toErange — 1 do k[1, i]: =Efun[i](x+h/2, inc (y, h*k[0, i]/2));

for i:=0 toErange — 1 do k[2, i]: =Efun[i](x+h/2, inc (y, h*k[1, i]/2));

for i:=0 toErange — 1 do k[3, i]: =Efun[i](x+h, inc (y, h*k[2, i]));

for i:=0 toErange — 1 do result[i]: =y[i]+h*(k[0, i] + 2* k[1, i]+ 2*k[2, i]+ k[3, i])/6;

end;

function TRunge. Value5(const x: real; const y: TRealArr; const h: real): TRealArr;

var i: integer; //счетчики.

k: array of TRealArr;

begin.

setlength (k, 6, Erange);

setLength (Result, Erange);

for i:=0 toErange — 1 dobegin.

k[0, i]: =EFun[i](x, y);

k[1, i]: =EFun[i](x + ¼*h, inc (y, h*k[0,i]/4));

k[2, i]: =EFun[i](x + 3/8*h, inc (y, h*(3/32*k[0,i] + 9/32*k[1,i])));

k[3, i]: =EFun[i](x + 12/13*h, inc (y, h*(1932/2197*k[0,i] - 7200/2197*k[1,i] + 7296/2197*k[2,i])));

k[4, i]: =EFun[i](x + h, inc (y, h*(439/216*k[0,i] - 8*k[1,i] + 3680/513*k[2,i] - 845/4104*k[3,i])));

k[5, i]: =EFun[i](x + h/2, inc (y, h*(- 8/27*k[0,i] + 2*k[1,i] - 3544/2565*k[2,i] + 1859/4104*k[3,i] - 11/40*k[4,i])));

result[i] :=y[i] + h*(16/135*k[0,i] + 6656/12 825*k[2,i] + 28 561/56430*k[3,i] +2/55*k[4,i]-9/50*k[5,i]);

end;

end;

//Проверка,.

functioncheckeps (consta, b: TRealArr): double;

var i: integer;

begin.

if length (a) length (b) then raiseException. Create ('Невернаядлинамассива. calc. CheckEps: a-'+inttostr (length (a))+', b-'+inttostr (length (b)));

result:=abs (a[0]-b[0]);

for i := 1 to high (a) do result:=max (result, abs (a[i] - b[i]));

end;

//Выполнение одного шага, со всеми проверками.

functionTRunge.calc (vlast, veps: double): TRes;

var.

h, hmin: double;

Arr: PResArr;

procedure Step;

vartmp, tmp2: TRes;

i: integer;

delta: double;

begin.

//Текущеезначениефункции.

tmp.arr:=value (Arr^[High (Arr^)]. x, Arr^[High (Arr^)].arr, h);

tmp.x:=Arr^[High (Arr^)]. x+h;

tmp2.arr:=value5(Arr^[High (Arr^)]. x, Arr^[High (Arr^)].arr, h);

tmp2.x:=tmp.x;

//Если функция по точности нам подходит.

delta:=checkeps (tmp2.arr, tmp. Arr)*abs (ELast-EFirst)/abs (hmin);

if delta.

hmin:=min (hmin, abs (h));

setlength (Arr^, length (Arr^)+1);

setlength (Arr^[High (Arr^)]. Arr, Erange);

for i := 0 to Erange-1 do.

Arr^[High (Arr^)]. Arr[i]:=tmp.arr[i];

Arr^[High (Arr^)]. x:=tmp2.x;

//Если точность хорошая, то увеличиваем длину шага.

if delta.

h:=2*h;

endelse //Если точность не достигнута, то уменьшаем шаг.

h:=h/2;

end;

begin.

Last:=vlast;

ifveps> 0 then.

Eeps:=veps;

h:=(Elast-Efirst)/10;

hmin:=abs (h);

ifElast>EFirstthenbegin.

Arr:=@Earr;

whileArr^[High (Arr^)]. x < last dobegin.

step;

end.

endelsebegin.

Arr:=@EarrNeg;

whileArr^[High (Arr^)]. x > last dobegin.

step;

end;

end;

result:=Arr^[High (Arr^)];

end;

дифференциальный аппроксимация массив численный В связи с некоторыми особенностями реализации и использованием сторонних компонентов построение фазового портрета в системе координат y (x) для поставленной задачи занимает продолжительное время.