Технический форум
Вернуться   Технический форум > Программирование > Форум программистов > Delphi, Kylix and Pascal


Ответ
 
Опции темы Опции просмотра
Старый 02.10.2019, 20:13   #1 (permalink)
t13ka
Новичок
 
Регистрация: 02.10.2019
Сообщений: 11
Сказал(а) спасибо: 0
Поблагодарили 0 раз(а) в 0 сообщениях
Репутация: 10
По умолчанию Найти решение задачи Коши для ОДУ методом Рунге-Кутта 4-го порядка

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

Старый 02.10.2019, 20:13
Helpmaster
Member
 
Аватар для Helpmaster
 
Регистрация: 08.03.2016
Сообщений: 0

Можно узнать множество полезной информации перейдя по этим ссылкам

Метод Рунге-Кутта 2-го порядка, Паскаль
Помогите найти интеграл и задачу коши операционным методом
Решение уравнения методом Рунге-Кутты 4-го порядка

Старый 03.10.2019, 20:15   #2 (permalink)
t13ka
Новичок
 
Регистрация: 02.10.2019
Сообщений: 11
Сказал(а) спасибо: 0
Поблагодарили 0 раз(а) в 0 сообщениях
Репутация: 10
По умолчанию

Я понял что в отличии от 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.
Миниатюры
5t9harjiorg.jpg  
t13ka вне форума   Ответить с цитированием
Старый 03.10.2019, 20:22   #3 (permalink)
t13ka
Новичок
 
Регистрация: 02.10.2019
Сообщений: 11
Сказал(а) спасибо: 0
Поблагодарили 0 раз(а) в 0 сообщениях
Репутация: 10
По умолчанию

ну или немного переделанный вариант, который работает, но выдает неправильные значения
Код:
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.
t13ka вне форума   Ответить с цитированием
Старый 03.10.2019, 20:40   #4 (permalink)
Vladimir_S
Специалист
 
Регистрация: 27.08.2008
Адрес: Санкт-Петербург
Сообщений: 27,807
Сказал(а) спасибо: 340
Поблагодарили 583 раз(а) в 208 сообщениях
Репутация: 113184
По умолчанию

Понятно.
То, что уравнений не два, а три — Вы поняли. Прекрасно. И даже тело программы записано правильно. Вот только... Понимаете, тут КАЖДАЯ из функций (у Вас — 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.
Миниатюры
fp01.jpg  
Vladimir_S вне форума   Ответить с цитированием
Старый 04.10.2019, 20:33   #5 (permalink)
t13ka
Новичок
 
Регистрация: 02.10.2019
Сообщений: 11
Сказал(а) спасибо: 0
Поблагодарили 0 раз(а) в 0 сообщениях
Репутация: 10
По умолчанию

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

Яндекс

Member
 
Регистрация: 31.10.2006
Сообщений: 40200
Записей в дневнике: 0
Сказал(а) спасибо: 0
Поблагодарили 0 раз(а) в 0 сообщениях
Репутация: 55070
Старый 04.10.2019, 22:23   #6 (permalink)
Vladimir_S
Специалист
 
Регистрация: 27.08.2008
Адрес: Санкт-Петербург
Сообщений: 27,807
Сказал(а) спасибо: 340
Поблагодарили 583 раз(а) в 208 сообщениях
Репутация: 113184
По умолчанию

Цитата:
Сообщение от t13ka Посмотреть сообщение
А как можно проверить правильность работы программы?
Даже и не знаю...
Цитата:
Сообщение от t13ka Посмотреть сообщение
Просто у меня проблемы с аналитическим решением этого уравнения
Ещё бы. Нелинейное третьего порядка!
Vladimir_S вне форума   Ответить с цитированием
Ads

Яндекс

Member
 
Регистрация: 31.10.2006
Сообщений: 40200
Записей в дневнике: 0
Сказал(а) спасибо: 0
Поблагодарили 0 раз(а) в 0 сообщениях
Репутация: 55070
Ответ

Опции темы
Опции просмотра

Ваши права в разделе
Вы не можете создавать новые темы
Вы не можете отвечать в темах
Вы не можете прикреплять вложения
Вы не можете редактировать свои сообщения

BB коды Вкл.
Смайлы Вкл.
[IMG] код Выкл.
HTML код Выкл.
Trackbacks are Вкл.
Pingbacks are Вкл.
Refbacks are Выкл.




Часовой пояс GMT +4, время: 22:10.

Powered by vBulletin® Version 6.2.5.
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.