Технический форум - Найти решение задачи Коши для ОДУ методом Рунге-Кутта 4-го порядка

Технический форум (http://www.tehnari.ru/)

- Delphi, Kylix and Pascal (http://www.tehnari.ru/f43/)

Найти решение задачи Коши для ОДУ методом Рунге-Кутта 4-го порядка

Найти решение задачи Коши для обыкновенного дифференциального уравнения методом Рунге-Кутты 4-го порядка точности.
Число разбиений n=20
y'''=x+x*y^2-(y')^2;
y(0)=1;
y'(0)=y''(0)=0
0<=x<= 2

Я понял что в отличии от y'' в которой у нас два уравнения, мое решение будут отличаться лишь добавлением еще одного уравнения. Тогда ниже приклепил свое решение на бумаге и код который я попытался переделать из вашей программы

Код:

Const

 x0 = 0;

 y0 = 1;

 z0 = 0;

 v0 = 0;

 h = 0.1;

 x_max = 2;



Var

 X, Y, Z, Y1, Z1: real;

 i,N:integer;



Function f(x, y, z: real): real;

begin

 f := x + x * y * y - z * z;

end;



Function g(x, y, z: real): real;

begin

 g := x + x * y * y - z * z;

end;



Function v(x, y, z: real): real;

begin

 v := z;

end;



Procedure RK(x, y, z: real; var Uy: real; var Uz: real; var Uv: real);

var k1, k2, k3, k4, q1, q2, q3, q4, p1, p2, p3, p4:real;

begin

 k1 := h * f(x, y, z);

 q1 := h * g(x, y, z);

 p1 := h * v(x, y, z);

 

 k2 := h * f(x + h / 2, y + q1 / 2, z + k1 / 2);

 q2 := h * g(x + h / 2, y + q1 / 2, z + k1 / 2);

 p2 := h * v(x + h / 2, y + q1 / 2, z + k1 / 2);

 

 k3 := h * f(x + h / 2, y + q2 / 2, z + k2 / 2);

 q3 := h * g(x + h / 2, y + q2 / 2, z + k2 / 2);

 p3 := h * v(x + h / 2, y + q2 / 2, z + k2 / 2);

 

 k4 := h * f(x + h, y + q3, z + q3);

 q4 := h * g(x + h, y + q3, z + q3);

 p4 := h * v(x + h, y + q3, z + q3);

 

 Uz := z + (k1 + 2 * k2 + 2 * k3 + k4) / 6;

 Uy := y + (q1 + 2 * q2 + 2 * q3 + q4) / 6;

 Uv := v + (p1 + 2 * p2 + 2 * p3 + p4) / 6;

end;



Begin

 N := Round((x_max - x0) / h);

 X := x0;

 Y := y0;

 Z := z0;

 V := v0;

 Writeln('x= ', X:4:2, '    y= ', Y:8:5);

 For i:=1 to N do

  begin

   RK(X, Y, Z, V, Y1 , Z1, V1);

   X := X + h;

   Y := Y1;

   Z := Z1;

   V := V1;

   Writeln('x= ', X:4:2, '    y= ', Y:8:5);

  end;

 Readln

End.

ну или немного переделанный вариант, который работает, но выдает неправильные значения

Код:

Const

 x0 = 0;

 y0 = 1;

 z0 = 0;

 v0 = 0;

 h = 0.1;

 x_max = 2;



Var

 X, Y, Z, Y1, Z1, C, C1: real;

 i,N:integer;



Function f(x, y, z: real): real;

begin

 f := x + x * y * y - z * z;

end;



Function g(x, y, z: real): real;

begin

 g := x + x * y * y - z * z;

end;



Function v(x, y, z: real): real;

begin

 v := z;

end;



Procedure RK(x, y, z: real; var Uy: real; var Uz: real; var Uv: real);

var k1, k2, k3, k4, q1, q2, q3, q4, p1, p2, p3, p4:real;

begin

 k1 := h * f(x, y, z);

 q1 := h * g(x, y, z);

 p1 := h * v(x, y, z);

 

 k2 := h * f(x + h / 2, y + q1 / 2, z + k1 / 2);

 q2 := h * g(x + h / 2, y + q1 / 2, z + k1 / 2);

 p2 := h * v(x + h / 2, y + q1 / 2, z + k1 / 2);

 

 k3 := h * f(x + h / 2, y + q2 / 2, z + k2 / 2);

 q3 := h * g(x + h / 2, y + q2 / 2, z + k2 / 2);

 p3 := h * v(x + h / 2, y + q2 / 2, z + k2 / 2);

 

 k4 := h * f(x + h, y + q3, z + q3);

 q4 := h * g(x + h, y + q3, z + q3);

 p4 := h * v(x + h, y + q3, z + q3);

 

 Uz := z + (k1 + 2 * k2 + 2 * k3 + k4) / 6;

 Uy := y + (q1 + 2 * q2 + 2 * q3 + q4) / 6;

 Uv := z + (p1 + 2 * p2 + 2 * p3 + p4) / 6;

end;



Begin

 N := Round((x_max - x0) / h);

 X := x0;

 Y := y0;

 Z := z0;

 C := v0;

 Writeln('x= ', X:4:2, '    y= ', Y:8:5);

 For i:=1 to N do

  begin

   RK(X, Y, Z, C, Y1 , Z1);

   X := X + h;

   Y := Y1;

   Z := Z1;

   C := C1;

   Writeln('x= ', X:4:2, '    y= ', Y:8:5);

  end;

 Readln

End.

Понятно.
То, что уравнений не два, а три — Вы поняли. Прекрасно. И даже тело программы записано правильно. Вот только... Понимаете, тут КАЖДАЯ из функций (у Вас — f, g и v) должна быть не от трёх, а от четырёх переменных. Кроме того, функция g записана неверно.
В общем, я тут нарисовал вариантик (немножко обозначения другие, в частности, у меня не v, а p) и решение задачки. Если не запутался в значках то, вроде, так. Попробуйте сравнить со своим кодом и разобраться:

Код:

Const

 x0=0.0;

 y0=1.0;

 z0=0.0;

 t0=0.0;

 N=20;

 x_max=2.0;



Var

 X,Y,Z,T,Y1,Z1,T1,h:real;

 i:integer;



Function f(x,y,z,t:real):real;

begin

 f:=x+x*y*y-z*z;

end;



Function g(x,y,z,t:real):real;

begin

 g:=t;

end;



Function p(x,y,z,t:real):real;

begin

 p:=z;

end;



Procedure RK(x,y,z,t:real; var Uy:real; var Uz:real; var Ut:real);

var k1,k2,k3,k4,q1,q2,q3,q4,r1,r2,r3,r4:real;

begin

 k1:=h*f(x,y,z,t);

 q1:=h*g(x,y,z,t);

 r1:=h*p(x,y,z,t);

 k2:=h*f(x+h/2,y+r1/2,z+q1/2,t+k1/2);

 q2:=h*g(x+h/2,y+r1/2,z+q1/2,t+k1/2);

 r2:=h*p(x+h/2,y+r1/2,z+q1/2,t+k1/2);

 k3:=h*f(x+h/2,y+r2/2,z+q2/2,t+k2/2);

 q3:=h*g(x+h/2,y+r2/2,z+q2/2,t+k2/2);

 r3:=h*p(x+h/2,y+r2/2,z+q2/2,t+k2/2);

 k4:=h*f(x+h,y+r3,z+q3,t+k3);

 q4:=h*g(x+h,y+r3,z+q3,t+k3);

 r4:=h*p(x+h,y+r3,z+q3,t+k3);

 Ut:=t+(k1+2*k2+2*k3+k4)/6;

 Uz:=z+(q1+2*q2+2*q3+q4)/6;

 Uy:=y+(r1+2*r2+2*r3+r4)/6;

end;



Begin

 h:=(x_max-x0)/N;

 X:=x0;

 Y:=y0;

 Z:=z0;

 T:=t0;

 Writeln('x = ',X:4:2,'    y = ',Y:8:5);

 For i:=1 to N do

  begin

   RK(X,Y,Z,T,Y1,Z1,T1);

   X:=X+h;

   Y:=Y1;

   Z:=Z1;

   T:=T1;

   Writeln('x = ',X:4:2,'    y = ',Y:8:5);

  end;

 Readln

End.

А как можно проверить правильность работы программы? Просто у меня проблемы с аналитическим решением этого уравнения

Цитата:

Сообщение от t13ka (Сообщение 2664399)

А как можно проверить правильность работы программы?

Даже и не знаю...

Цитата:

Сообщение от t13ka (Сообщение 2664399)

Просто у меня проблемы с аналитическим решением этого уравнения

Ещё бы. Нелинейное третьего порядка!