Погрешность метода Рунге-Кутта

← →
ksevelyar (2009-04-15 19:09) [0]

Вроде родил алгоритм

program runge_kutta; var x, y, h: real; function f(x,y: real): real; begin f := x - y*y end; function rk(x, y, h: real): real; var m: array [0..3] of real; h_2: real; begin h_2 := 0.5 * h; m[0] := f(x,y); m[1] := f(x + h_2, y + h_2*m[0]); m[2] := f(x + h_2, y + h_2*m[1]); m[3] := f(x + h, y + h*m[2]); rk := y + h*(m[0] + 2*m[1] + 2*m[2] + m[3])/6; end; begin x := 0; y := 0; h := 0.01; writeln("x":11, "y":11); while x <= 0.1 do begin writeln(x:11:6, y:11:6); y := rk(x, y, h); x := x + h; end; end.

Но намертво застрял с вычислением погрешности - например eps = 0,0001

Буду очень признателен любой информации.

← →
Denis__ © (2009-04-15 20:14) [1]

В описании алгоритма должно быть написано(А то что Вы написали - это не алгоритм. Это программа на Паскале, возможно реализующая выполнение алгоритма). И потом Вы не указали, метод вычисления чего?

← →
Юрий Зотов © (2009-04-15 22:55) [2]

> ksevelyar (15.04.09 19:09)

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

> Denis__ © (15.04.09 20:14) [1]

> то что Вы написали - это не алгоритм. Это программа на Паскале
Программа на любом языке - это тоже один из способов описания алгоритма.

> Вы не указали, метод вычисления чего?
Вы в курсе, что такое метод Рунге-Кутта? Если да, то вопросов возникать не должно, а если нет - то возникают другие вопросы.

← →
Юрий Зотов © (2009-04-15 23:33) [3]

> ksevelyar (15.04.09 19:09)

Вдогонку: еще можно использовать метод Адамса (предсказания-коррекции), там погрешность оценивается непосредственно.

← →
Григорьев Антон © (2009-04-16 10:29) [4]

Есть два способа оценки погрешнгости методов Рунге-Кутта: правило Рунге и контрольный член Егорова. Гуглите по этим словам.

← →
Jeer © (2009-04-16 10:37) [5]

> ksevelyar (15.04.09 19:09)
>
> Вроде родил алгоритм

Для начала - сей модуль родил не ты, а отсюда и непонятки у тебя :)

This archive contains the following files. They are all concerned
with some aspect of numerical analysis. In most cases the
programs include a reference to their source. Perhaps the most
generally useful and understandable text is "Applied Numerical
Analysis", Curtis F. Gerald -- his sample programs are about as
structured as it is possible to be in Fortran. Harley Flanders"
"Scientific Pascal" has excellent sample programs, but is not
very strong on the theory behind them. Whether you feel like
looking up the theory or not, it usually a good idea to test the
program by running it on a case where the analytic answer is
known.

LIKKER.PAS Solves a heat transfer problem
LINSYS.PAS System of linear equations
MATLIB.PAS Collection matrix operations. Should be used with
discretion - calculating the inverse is NOT a good
way to solve a linear system.
MULTINT.PAS Multiple definite integration
NEWTON.PAS Finds a zero of a function
PANREIF2.PAS Invert a real matrix
RK.PAS Runge-Kutta solution of ordinary differential
equation. (o.d.e.)
RKF45.PAS Adaptive Runge-Kutta-Fehlberg (~5th order accuracy)
RKNRLC.PAS Runge-Kutta-Nystrom -- 2nd. order equations
RKSYS.PAS System of o.d.e."s
RKV71.PAS Runge-Kutta-Verner, o.d.e"s with 7th order accuracy
...

program runge_kutta;
var x, y, h: real;

function f(x,y: real): real;
begin f := x - y*y end;

function rk(x, y, h: real): real;
var m: array [0..3] of real;
h_2: real;
begin
h_2 := 0.5 * h;
m[0] := f(x,y);
m[1] := f(x + h_2, y + h_2*m[0]);
m[2] := f(x + h_2, y + h_2*m[1]);
m[3] := f(x + h, y + h*m[2]);
rk := y + h*(m[0] + 2*m[1] + 2*m[2] + m[3])/6;
end;

begin
x := 0; y := 0; h := 0.01;
writeln(" Runge-Kutta Method");
writeln(" y"" = x - y^2, y(0) = 0");
writeln("x":11, "y":11);
while x <= 0.1 do begin
writeln(x:11:6, y:11:6);
y := rk(x, y, h);
x := x + h;
end;
end.

Submitted by MARK STEPHEN.
Leave a message if you have any problems with this stuff, but any
numerical analysis text will probably be more helpful than I can.

********************

Это метод Р-К четвертого порядка
Оценку погрешности можно производить методом Рунге ( двойное вычисление с разным шагом )

Для указанного метода при кратности шага k=2h имеем оценку погрешности

Ro = [ Yh(x) - Yk(x)]/15

Для оценки погрешности метода Р-К 4-го порядка на каждом шаге существует его разновидность, предложенная Мерсеном.
Метод допускает автоматизацию выбора шага интегрирования.

← →
ksevelyar (2009-04-16 14:26) [6]

Спасибо всем кто ответил! Да, код написал не я - я лишь адаптировал его к дельфи (сюда выложил именно вариант паскаля с умолчательной функцией чтобы не отвлекать внимание).

>
> Для указанного метода при кратности шага k=2h имеем оценку
> погрешности
>
> Ro = [ Yh(x) - Yk(x)]/15
>
> Для оценки погрешности метода Р-К 4-го порядка на каждом
> шаге существует его разновидность, предложенная Мерсеном.
>
> Метод допускает автоматизацию выбора шага интегрирования.
>

То есть как я понял:

1) считается значение у при h, считается значение y при h2:= h/2
2) проверяется погрешность Ro = [Yh2(x) - Yh(x)]/15 (тут модуль или просто скобки?)
Если Ro <= eps то y:=yh(x), иначе h:=h/2; и к пункту 1

Исправил код на следующий, и всё равно не работает, массивы забиваются нулями и всё :(

eps:=0.00001; x := 0; t[0]:=0; y[0] := 0; h := 0.001; n:=1; repeat yh1 := rk(x, yh1, h); yh2 := rk(x, yh2, h/2); Ro:= abs(Yh2 - Yh1)/15; if Ro < eps then begin y[n]:=yh1; t[n]:=x; n:=n+1; x:=x+h end else h:=h/2; until x <= 0.1;

← →
ksevelyar (2009-04-16 14:31) [7]

В алгоритме похоже ошибка: надо until x >= 0.1;

нов ведь и в оригинале было while x <= 0.1 do begin...

← →
ksevelyar © (2009-04-16 14:55) [8]

И тем не менее - при каждом нажатии на кнопку вычислить результат увеличивается :(

Массивы очищать пробовал:

> procedure TForm1.Button1Click(Sender: TObject);
>
> begin
> for n := 1 to 100 do
> begin
> y[n]:=0; t[n]:=0;
> end;
>
>
> eps:=0.00001; Ro:=0;
> x := 0; h := 0.01; n:=1;
> repeat
> yh1 := rk(x, yh1, h);
> yh2 := rk(x, yh2, h/2);
> Ro:= (Yh2 - Yh1)/15;
> if Ro < eps then
> begin
> y[n]:=yh1;
>
> t[n]:=x;
>
> n:=n+1;
>
> x:=x+h
> end
> else h:=h/2;
>
> until x >= 0.1;
> for n:= 1 to 100 do begin
> StringGrid1.Cells[1,n]:= FloatToStr(y[n]);
> StringGrid1.Cells[0,n]:= FloatToStr(t[n]);
> end;

Надо было обнулять ещё промежуточные у:
Ro:=0; yh1:=0; yh2:=0; x := 0; h := 0.01; n:=0; while x <= 0.1 do begin if Ro < eIL1 then begin y[n]:=yh1; t[n]:=x; n:=n+1; x:=x+h end else h:=h/2; yh1 := rk(x, yh1, h); yh2 := rk(x, yh2, h/2); Ro:= (Yh2 - Yh1)/15; end;

Однако с разной погрешностью результат не меняется вообще. Значит я что-то сделал не так...

Извиняюсь за то, что создааю новые сообщения - не могу понять как править старые (я зарегестрировался).

> не могу понять как править старые

Никак. Такой возможности нет.

Погрешность метода Рунге-Кутта Найти похожие ветки