Метод Рунге кутта для системы из 2х уравнений Mathematica

[Alexeybk5]

Здравствуйте, уважаемые товарищи.

Подскажите пожалуйста, как можно с помощью Mathematica реализовать метод рунге кутта для системы из 2 уравнений с двумя переменными

$\frac{dx}{dt}=ax-bxy$
$\frac{dy}{dt}=-cy+dxy$

где мы знаем значения $a=2,\,\,\,\,\,b=0.002,\,\,\,\,\,\,c=1,\,\,\,\,\,\,\, d=0.002$

В итоге надо применить метод к системе в таком ввиде


$\frac{dx}{dt}=2x-0,002xy$
$\frac{dy}{dt}=-1y+0,002xy$


Я нашёл решение похожей задачи и код для неё тут
http://wps.prenhall.com/wps/media/objects/884/905485/chapt4/proj4.3A/proj4-3A.pdf
В коде реализован метод Рунге кутта для системы из 2 уравенений, для примера
$x(t)=\cos(\pi t),\,\,\,\,\,\,\y(t)=\sin(\pi t)$
$x'=-ay,\,\,\,\,\,\,\,\,\, y' = bx$

pi = N[Pi]; 

f[t_,x_,y_] := -pi*y 
g[t_,x_,y_] :=  pi*x 


t = t0; x = x0; y = y0;

Do[k1 = f[t, x, y];
 l1 = g[t, x, y];(*left-hand slopes*)
 k2 = f[t + h/2, x + h*k1/2, y + h*l1/2];
 l2 = g[t + h/2, x + h*k1/2, y + h*l1/2];(*1st midpt slopes*)
 k3 = f[t + h/2, x + h*k2/2, y + h*l2/2];
 l3 = g[t + h/2, x + h*k2/2, y + h*l2/2];(*2nd midpt slopes*)
 k4 = f[t + h, x + h*k3, y + h*l3];
 l4 = g[t + h, x + h*k3, y + h*l3];(*right-hand slopes*)
 k = (k1 + 2*k2 + 2*k3 + k4)/6;(*average x-slope*)
 l = (l1 + 2*l2 + 2*l3 + l4)/6;(*average y-slope*)
 x = x + h*k;(*update x*)y = y + h*l;(*update y*)t = t + h;(*update t*)
 If[ Floor[i/m] == i/m, 
       Print[180*t,PaddedForm[x,{10,4}], 
       PaddedForm[y,{10,4}]]], 
(* display each mth value *) 
   {i,1,n} ] 


15.      0.9659      0.2588 
30.      0.8660      0.5000 
45.      0.7071      0.7071 
60.      0.5000      0.8660 
75.      0.2588      0.9659 
90.      0.0000      1.0000 

Я попытался переделать, но в результате программа выдаёт очень странный результат

f[t_, x_, y_] := 2*x - 0.002*x*y
g[t_, x_, y_] := -1*y + 0.002*x*y

t0 = 0; x0 = 1; y0 = 0; tf = 0.5;

n = 12;(*number of subintervals*)m = 2;
(*to print every mth step*)h = (tf - t0)/n;(*step size*)

t = t0; x = x0; y = y0;

Do[k1 = f[t, x, y];
 l1 = g[t, x, y];(*left-hand slopes*)
 k2 = f[t + h/2, x + h*k1/2, y + h*l1/2];
 l2 = g[t + h/2, x + h*k1/2, y + h*l1/2];(*1st midpt slopes*)
 k3 = f[t + h/2, x + h*k2/2, y + h*l2/2];
 l3 = g[t + h/2, x + h*k2/2, y + h*l2/2];(*2nd midpt slopes*)
 k4 = f[t + h, x + h*k3, y + h*l3];
 l4 = g[t + h, x + h*k3, y + h*l3];(*right-hand slopes*)
 k = (k1 + 2*k2 + 2*k3 + k4)/6;(*average x-slope*)
 l = (l1 + 2*l2 + 2*l3 + l4)/6;(*average y-slope*)
 x = x + h*k;(*update x*)y = y + h*l;(*update y*)t = t + h;(*update t*)
 If[Floor[i/m] == i/m, 
  Print[180*t, PaddedForm[x, {10, 4}], 
   PaddedForm[y, {10, 4}]]],(*display each mth value*){i, 1, n}]

(Оффтоп)

15. 0.9659 0.2588
30. 0.8660 0.5000
45. 0.7071 0.7071
60. 0.5000 0.8660
75. 0.2588 0.9659
90. 3.8197*10^(-6) 1.0000

Помогите пожалуйста понять, что нужно изменить
Спасибо большое

[Yu_K]

А в чем странность?

[Alexeybk5]

Сложность в том, что я не пойму, какой период выбрать ? Это к математической модели "хищник-жертва"

[Yu_K]

Поделите одно на другое и посмотрите траектории в плоскости $(x,y)$

$$\frac{dx}{dy}=\frac{2x-0,002xy}{-1y+0,002xy}$$

- будут там замкнутые траектории? Вот если поменять местами правые части Ваших двух уравнений - то похоже будут замкнутые трактории - ну или хотя бы фокусы - а тут похоже нет их.