I'm trying to approximate a systems of differential equations using Euler's method, but at some point one of the variables is evaluated to NaN and from that point they all become Nan's. I've tried debugging and the onset of this varies depending on the choice of parameters (especially initial conditions and step size). For the below code z = NaN at the 989th iteration.
I tried extracting the values of the variables at that point and calculating the next step separately and it turns out they work fine. Actually, I can continue evaluating the functions for another 500ish steps until the same thing happens.
h = 0.05; %step
final = 100; %simulation time
t = 0:h:final;n = 1; %Hill coefs
m = 1;p = 1;v0 = 0.000001; v1 = 0.000003; B = 0.01; VM2 = 0.0001; VM3 = 0.001;K2 = 0.000001;Kr = 0.0001;Ka = 0.0000025;kf = 0.1;k = 2;y = 24; %variable initial value
z = 0.1; %variable initial valuefor n = 1:final/hz(n+1) = z(n) + h*(v0 + v1*B - VM2*ratio(n,z(n),K2) + VM3*ratio(m,y(n),Kr)*ratio(p,z(n),Ka) + kf*y(n) - k*z(n));y(n+1) = y(n) + h * (VM2*ratio(n,z(n),K2) - VM3*ratio(m,y(n),Kr)*ratio(p,z(n),Ka)- kf*y(n));end
The function ratio() comes from a separate file, so I add it here:
function [ factor ] = ratio( ex, var, K )factor = power(var,ex) / ( power(K,ex) + power(var,ex) );end
I expect this system to approach a stable limit cycle, so both y and z should keep oscillating forever. It might turn out that the parameters not chosen appropriately, but the NaN problem is the main issue now.
Since I'm a beginner, any additional comments on my "bad programming practices" are welcome.
Best Answer