Метод Рунге-Кутта-Мерсона для дифференциальных уравнений 2 порядка

@Taras_lebdyk · Регистрация: 17.09.2013

Студворк — интернет-сервис помощи студентам

Добрый день. Помогите пожалуйста понять в чем ошибка. У меня дано такое уравнение:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?0.3{y}^{''}=\frac{0.15y}{x}-\frac{{y}^{'}}{{x}^{2}}$
Начальные условия:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?y(1.0)=0; {y}^{'}=1.0; n=60; a=1; b=5;$
У меня есть программа для дифференциальных уравнений 1 порядка:

Pascal

Program vvv; 
const
e=0.0001;
var x0,xn,y0,h,x,y,k5,k1,k2,k3,k4,a,b,r:real; 
    n,i:integer; 
 
function F(x,y:real):real; 
begin 
   F:=0.25*y*y*(exp(x-1)-exp(-x+1))/(exp(x-1)+exp(-x+1))-x*x*y;  
end; 
begin 
a:=0.36;
b:=1.8;
n:=150;
x0:=0.36;
y0:=1.5;
i:=0;
h:=(b-a)/n; 
x:=x0; y:=y0; 
repeat 
     k1:=h*F(x,y);  
     k2:=h*F(x+h/3,y+k1/3); 
     k3:=h*F(x+h/3,y+k1/6+k2/6);  
     k4:=h*F(x+h/2,y+k1/8+3*k3/8);
     k5:=h*F(x+h,y+k1/2-3*k3/2+2*k4);
     r:=(-2*k1-9*k3-8*k4+k5)*h/30;
     if abs(r)>e then h:=h/2;
     if abs(r)<e/30 then h:=h*2;
     x:=x+h;
     y:=y+(k1+4*k4+k5)/6;
     i:=i+1;
        
     Writeln('|', i:2, ' |', x:5:3, ' |', y:7:3, ' |');
 until i>n;
readln; 
end.

Я переделай для 2 порядка но у меня не правильно работает( незнаю в чем ошибка подскажите пожалуйста. Вот код программы:

Pascal

Program vvv; 
var x0,xn,y0,h,x,y,k1,k2,k3,k4,k5,q1,q2,q3,q4,q5,r,z,z0,a,b:real; 
    n,i:integer; 
const
e=0.0001;
function F(x,y,z:real):real; 
begin 
   F:=z; 
end; 
function G(x,y,z:real):real; 
begin 
   G:=(((0.15*y)/(x*0.3))-(z/((x*x)*0.3))); 
end; 
begin 
x0:=1;
y0:=0;
z0:=1;
a:=1;
b:=5;
n:=60;
h:=(b-a)/n; 
x:=x0; y:=y0; z:=z0;
i:=0;
repeat    
     k1:=h*F(x,y,z); 
     q1:=h*G(x,y,z); 
     k2:=h*F(x+h/3,y+k1/3,z+q1/3);
     q2:=h*G(x+h/2,y+k1/3,z+q1/3); 
     k3:=h*F(x+h/3,y+k1/6+k2/6,z+q1/6+q2/6);  
     q3:=h*G(x+h/3,y+k1/6+k2/6,z+q1/6+q2/6); 
     k4:=h*F(x+h/2,y+k1/8+3*k3/8,z+q1/8+3*q3/8); 
     q4:=h*G(x+h/2,y+k1/8+3*k3/8,z+q1/8+3*q3/8);
     k5:=h*F(x+h,y+k1/2-3*k3/2+2*k4,z+q1/2-3*q3/2+2*q4);
     q5:=h*G(x+h,y+k1/2-3*k3/2+2*k4,z+q1/2-3*q3/2+2*q4);
      r:=(-2*k1-9*k3-8*k4+k5)*h/30;
     if abs(r)>e then h:=h/2;
     if abs(r)<e/30 then h:=h*2;
     x:=x+h;
     y:=y+(k1+4*k4+k5)/6;
     z:=z+(q1+4*q4+q5)/6;
     i:=i+1;
    
     
     Writeln('|', i:2, ' |', x:5:3, ' |', y:7:3, ' |', z:7:3, ' |');
until i>n; 
readln; 
end.

ФедосеевПавел · 30.03.2015, 23:10

Не по теме:

Странная ситуация. Поиск в интернете даёт по разным ссылкам разные варианты формул метода и оценки погрешности на текущем шаге. Я смог запустить, описанный здесь. Некоторые другие давали расходящиеся результаты для теста y''=-sin(x).

А по вопросу. Если всё же почитать описание метода, то можно заметить, что это метод с адаптацией шага интегрирования (строки 35-37). И значит, что условие завершения должно быть не i>n, а x>b. А условие печати строки x>i*(b-a)/n, а также вместе с печатью увеличивать i на 1. Может быть ещё добавить условие, чтобы шаг не увеличивался больше (b-a)/n.

ФедосеевПавел Модератор 8317 / 4216 / 1602 Регистрация: 01.02.2015 Сообщений: 13,132 Записей в блоге: 4
	30.03.2015, 23:10	2
	Не по теме: Странная ситуация. Поиск в интернете даёт по разным ссылкам разные варианты формул метода и оценки погрешности на текущем шаге. Я смог запустить, описанный здесь. Некоторые другие давали расходящиеся результаты для теста y''=-sin(x). А по вопросу. Если всё же почитать описание метода, то можно заметить, что это метод с адаптацией шага интегрирования (строки 35-37). И значит, что условие завершения должно быть не i>n, а x>b. А условие печати строки x>i*(b-a)/n, а также вместе с печатью увеличивать i на 1. Может быть ещё добавить условие, чтобы шаг не увеличивался больше (b-a)/n. 0

Опции темы