Ode-system output gives null values in Matlab

Question

I am trying to compute the numerical solution (over the interval [0,2]) for this system of differential equations using the Runga Katta (4th Order) Method and a step of my choice.

% This is not just equations, not MATLAB code
dy1/dt = 1.3 * (y2-y1) + 1.04 * 10^4 * k * y2
dy2/dt = 1.88 * 10^3 * (y4-(1+k))
dy3/dt = 1752+(266.1 * y2)-(269.3 * y3)
dy4/dt = 0.1+(320 * y2)-(321 * y4)
k = 6 * 10^(-4) * exp(20.7 - 15 * 10^3/y1)
y0 = transpose([759.167,0,600,0.1])

I adopted the following code, from this answer which ran smoothly but output a table of zeroes (barring a lone entry). I tried different step-sizes but to no avail.

clc

h = 0.2; % Step size.
t = (0:h:2)'; 
N = length(t);

y1 = zeros(N,1);    
y2 = zeros(N,1);    
y3 = zeros(N,1);
y4 = zeros(N,1);

% Initial conditions
y1(1) = 759.167;     
y2(1) = 0;    
y3(1) = 600;
y4(1) = 0.1;

% Initialise derivative functions
dy1 = @(t, y1, y2, y3, y4) (1.3*(y2-y1))+1.04*10^4*6*10^(-4)*exp(20.7 - (15*10^(3)/y1))*y2;
dy2 = @(t, y1, y2, y3, y4) 1.88*10^3*(y4-(1 + 6*10^(-4)*exp(20.7 - (15*10^(3)/y1))));
dy3 = @(t, y1, y2, y3, y4) 1752+(266.1*y2)-(269.3*y3);
dy4 = @(t, y1, y2, y3, y4) 0.1+(320*y2)-(321*y4);

% Initialise K vectors
ky1 = zeros(1,4);
ky2 = zeros(1,4);
ky3 = zeros(1,4);
ky4 = zeros(1,4);

% RK4 coefficients
b = [1 2 2 1];

% Iterate, computing each K value in turn, then the i+1 step values
for i = 1:(N-1)        
    ky1(1) = dy1(t(i), y1(i), y2(i), y3(i), y4(i));
    ky2(1) = dy2(t(i), y1(i), y2(i), y3(i), y4(i));
    ky3(1) = dy3(t(i), y1(i), y2(i), y3(i), y4(i));
    ky4(1) = dy4(t(i), y1(i), y2(i), y3(i), y4(i));

    ky1(2) = dy1(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));
    ky2(2) = dy2(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));
    ky3(2) = dy3(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));
    ky4(2) = dy4(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));

    ky1(3) = dy1(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));
    ky2(3) = dy2(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));
    ky3(3) = dy3(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));
    ky4(3) = dy4(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));

    ky1(4) = dy1(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));
    ky2(4) = dy2(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));
    ky3(4) = dy3(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));
    ky4(4) = dy4(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));

    %y1(i+1)
    y1(i+1) = y1(i) + (h/6)*sum(b.*ky1);       
    y2(i+1) = y2(i) + (h/6)*sum(b.*ky2);       
    y3(i+1) = y3(i) + (h/6)*sum(b.*ky3);
    y4(i+1) = y4(i) + (h/6)*sum(b.*ky4);
end

txyz = [t,y1,y2,y3]

Please post your problem as text in your question and highlight your *specific* issue, at this minute this question is 90% just a copy/paste of my code from that other answer with little context on what you've changed or why... It would also probably be easier if you used the standalone RK4 function which I included in that answer at the bottom, rather than adapting my adaptation of that question's code! — Wolfie, Jan 06 '21 at 14:52
This stack-exchange does not support Latex. and I am not able write the question in a better manner. My question simply adds another variable (and another equation). I have accordingly scaled up your first code. My specific issue is the code does not flag any error but returns a table filled with zeroes. — Alpha Range Thirty, Jan 06 '21 at 15:22
I followed your standalone function as well but it shew me "Not Enough Input Arguments". Code over [here](https://pastebin.com/q0GzYFme) — Alpha Range Thirty, Jan 06 '21 at 15:53
Your equations seem to be completely independent of `t`, how do you expect them to evolve over time (which the Runga Kutta method aims to give you here)? That seems like an issue, and then running these equations anyway with the generic `RK4` function, you can see that the `y` values blow up to `1e+299` before becoming `NaN` because the values are too large to store at double precision. You can see from some back-of-the-envelope maths that the `y` values will blow up quickly from their non-linear dependencies and large coefficients. — Wolfie, Jan 06 '21 at 17:49
Why do we have a problem, even if the _derivatives_ (**not** variables) are completely independent of `t`? — Alpha Range Thirty, Jan 11 '21 at 12:57
In this case, the variables blow up, as I mentioned, I thought maybe you were missing a `t` dependency which would prevent that. Do you see what I mean about it being clear your variables will blow up? — Wolfie, Jan 11 '21 at 17:17

Ode-system output gives null values in Matlab

0 Answers0