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

Размещено на http://www.allbest.ru

1. Теория

Методы Рунге-Кутты являются модифицированными методами Эйлера. Они служат для численного решения обыкновенных дифференциальных уравнений и их систем. Наиболее часто реализуется и используется метод Рунге-Кутты 4 порядка точности. Построение схем более высокого порядка приводит к громоздким расчетам и сопряжено с большими вычислительными трудностями. Схема Эйлера, является представителем семейства схем Р-К. Построение схем основано на разложении функции в ряд Тейлора вблизи инвариантной точки, и взятие конечного числа членов разложения.

Пусть дана задача:

.

Тогда приближенное решение задается формулой:

,

где h - величина шага сетки, а kвычисляется по формуле:

.

.

.

Метод определяется числомs и коэффициентами b_i,a_i,jи c_i

Причем:

.

Порядок аппроксимации (точнее порядок сходимости) определить сложно, однако можно определить количество этапов, необходимых для реализации метода.

.

Конкретно метод Р-К 4 для системы ДУ

.

имеет следующий вид.

.

.

.

.

.

По сути, берется значение методом Эйлера, и уточняется 3 раза. За счет этого возрастает точность вычислений, ценой вычислительных трудностей.

Для метода Рунге-Кутты 5 порядка уравнения будут такими:

.

.

.

.

.

.

.

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

Метод Рунге-Кутты-Фельберга, заключается в том, что на каждом шаге метода точность функции определяется разностью значений между результатами методов РК-4 и РК-5 (Поэтому этот метод иногда называют РК-45). Если они отличаются не более чем на е - локальную погрешность, то значение, уточненное по Рунге, считается приближенным значением функции в точке на рассматриваемом шаге.

2. Реализация

Программа, реализующая адаптивную процедуру Рунге-Кутты 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 + 1/4*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/12825*k[2,i] + 28561/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 <Eepsthenbegin

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 <Eeps/1000then

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) для поставленной задачи занимает продолжительное время.

Размещено на Allbest.ru

...