Решение ОДУ 2-г порядка методом Рунге-Кутты 4-го порядка

[Imaginarium]

Здравствуйте.

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

Задача: есть система частиц (тел, заряженных шариков), которые располагаются вдоль горизонтальной линии (координата $x$) никого уровня как бусины через равные промежутки. Двигаться по горизонтали тела не могут, только вверх-вниз (по координате $r$) . На каждое тело действуют все стоящие справа от него тела (соответственно, самое правое все время остается неподвижным) по некоему закону, определяя его положение: $\frac {d^2r^{i}} {dt^2} = \sum _{j=1}^{i-1} A \cdot r^i(t)\cdot \sin [k(x_i-x_j)]$, где $i$- номер тела;
с начальными условиями $\begin{cases} r(0) = r_0\\ r'(0) = 0\\ \end{cases}\\ $.
Для того, чтобы применить метод, я делаю из ОДУ 2-го порядка систему 1-го порядка:
$\begin{cases} y' = y_1\\ y'_1 = Cy\\ \end{cases}\\ C=\operatorname{const} вычисляется для каждого тела заново, каждый шаг по времени $.
Теперь надо написать выражения для всех коэффициентов метода, и т.к. у нас 4-й порядок точности и 2 уравнения, то коэффициентов будет 4 штуки для каждого уравнения (k и m).
И вот тут-то я и начинаю путаться-
выписываем первую пару коэффициентов:
$k_1 = h\cdot 0 , m_1 = h \cdot r_0$
А вторая пара? К чему мне прибавлять половину шага?

Хочу написать программу на Matlab, не пользуясь встроенными функциями типа ode45.

Спасибо.

[Bauyrzhan]

попробуй это, только не забудь подставить свою функцию))

clear,clc;
close all;
h=0.1;
x=0:h:4;
z=zeros(1,length(x));
y=zeros(1,length(x));
F_xy=@(z,x,y)z;
G_xy=@(x,y,z)1+exp(x)-y;  %функция
y(1)=2.5;                 %начальное условие
z(1)=1.5;                 %нач. усл.  
n=length(x)-1;
for i=1:n
    K1=F_xy(x(i),y(i),z(i));       
    q1=G_xy(x(i),y(i),z(i));
    K2=F_xy(x(i)+h/2,y(i)+h/2*K1,z(i)+h/2*q1);
    q2=G_xy(x(i)+h/2,y(i)+h/2*K1,z(i)+h/2*q1);
    K3=F_xy(x(i)+h/2,y(i)+h/2*K2,z(i)+h/2*q2);
    q3=G_xy(x(i)+h/2,y(i)+h/2*K2,z(i)+h/2*q2);
    K4=F_xy(x(i)+h,y(i)+h*K3,z(i)+h*q3);
    q4=G_xy(x(i)+h,y(i)+h*K3,z(i)+h*q3);   
    z(i+1)=z(i)+(h/6)*(K1+2*K2+2*K3+K4);   
    y(i+1)=y(i)+(h/6)*(q1+2*q2+2*q3+q4);
end
plot(x,y,x,z,'Linewidth',2);grid;
%figure(1); plot(x,y,'r');
%figure(2); plot(x,z,'g');
legend('y''','y''''',2);

[Deggial]

Bauyrzhan в сообщении #722582 писал(а):

попробуй это

Bauyrzhan, замечание за фамильярность, на форуме следует обращаться друг к другу на "Вы".

Правила форума в сообщении #722582 писал(а):

1) Нарушением считается:
...
е) Провокационные и вызывающие сообщения, фамильярность (у нас принято обращаться друг к другу на "Вы")...