Не работает реализация Рунге-Кутты 4 порядка MATHLAB

Question

В первой итерации рассчитываются шаги k1,k2,k3,k4, но как только я добавляю последний шаг формулы метода РК X(i+1) = X(i) + 1/6.*(k1X + 2.*k2X + 2.*k3X + k4X) все ломается и в значениях k появляется NaN. Подозреваю, что проблема в переполнении, но как правильно определить функции не понимаю. Вывод, конечно, будет для всех функций, но пока сделала только для Х чтобы легче было ориентироваться по графику. Вот код на MATHLAB:

 clc;                                               % Clears the screen
clear; 

%coefficients
v = 0.4;
kx = 0.4;
ky = 0.4;
kz = 0.4;
kg = 0.4;
kf = 0.5;
ax = 0.4;
ay = 0.4;
agf = 0.06;
hz = 0.9;
m0 = 1;


%step size
h = 0.04;

tfinal = 100;
n = ceil(tfinal/h);
t = 0:h:100;



%initail conditions
X(1) = 1.0;
Y(1) = 0.0;
Z(1) = 0.0;
G(1) = 5.0;
F(1) = 0.0;



%define function 
fx = @(t,X,G) v-ax*X*G-kx*X;

fy = @(t,X,Y,G) ax*X*G+m0*Y-ay*Y*G- ky*Y;

fz = @(t,Y,G,Z) ay*Y*G-kz*Z;

fg = @(t,F,G) -kg*G-agf*G*F;

ff = @(t,F,G,Z) hz*Z-agf*G*F- kf*F;

 

   

 for i = 1:n 

   
    %KuttaSteps
     
    k1X = h*fx(t(i), X(i), G(i));
    k1Y = h*fy(t(i), X(i), Y(i), G(i));
    k1Z = h*fz(t(i), Y(i), G(i), Z(i));
    k1G = h*fg(t(i), F(i), G(i));
    k1F = h*ff(t(i), F(i), G(i), Z(i));


    k2X = h*fx(t(i)+h/2., X(i)+k1X/2., G(i)+k1G/2.);
    k2Y = h*fy(t(i)+h/2., X(i)+k1X/2., Y(i)+k1Y/2., G(i)+k1G/2.);
    k2Z = h*fz(t(i)+h/2., Y(i)+k1Y/2., G(i)+k1G/2., Z(i)+k1Z/2.);
    k2G = h*fg(t(i)+h/2., F(i)+k1F/2., G(i)+k1G/2.);
    k2F = h*ff(t(i)+h/2., F(i)+k1F/2., G(i)+k1G/2., Z(i)+k1Z/2.);


    k3X = h*fx(t(i)+h/2., X(i)+k1X/2., Y(i)+k2Y/2.);
    k3Y = h*fy(t(i)+h/2., X(i)+k1X/2., Y(i)+k2Y/2., G(i)+k2G/2.);
    k3Z = h*fz(t(i)+h/2., Y(i)+k2Y/2., G(i)+k2G/2., Z(i)+k2Z/2.);
    k3G = h*fg(t(i)+h/2., F(i)+k2F/2., G(i)+k2G/2.);
    k3F = h*ff(t(i)+h/2., F(i)+k2F/2., G(i)+k2G/2., Z(i)+k2Z/2.);

    k4X = h*fx(t(i)+h, X(i)+k1X, Y(i)+k3Y);
    k4Y = h*fy(t(i)+h, X(i)+k1X, Y(i)+k3Y, G(i)+k3G);
    k4Z = h*fz(t(i)+h, Y(i)+k3Y, G(i)+k3G, Z(i)+k3Z);
    k4G = h*fg(t(i)+h, F(i)+k3F, G(i)+k3G);
    k4F = h*ff(t(i)+h, F(i)+k3F, G(i)+k3G, Z(i)+k3Z);
    
    X(i+1) = X(i) + 1/6.*(k1X + 2.*k2X + 2.*k3X + k4X);
    Y(i+1) = Y(i) + 1/6.*(k1Y + 2.*k2Y + 2.*k3Y + k4Y);
    Z(i+1) = Z(i) + 1/6.*(k1Z + 2.*k2Z + 2.*k3Z + k4Z);
    G(i+1) = G(i) + 1/6.*(k1G + 2.*k2G + 2.*k3G + k4G);
    F(i+1) = F(i) + 1/6.*(k1F + 2.*k2F + 2.*k3F + k4F);
     
     
        
 end
 
 
 %plotting
 figure(1); clf(1)
 plot(t,X)
 xlabel('Time')
 ylabel('X')

БЛОГ НА HUSL

Не работает реализация Рунге-Кутты 4 порядка MATHLAB

Ответы (0 шт):