For loops breaks into NaN's

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 value

for n = 1:final/h
z(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.

Answer 1

[David Goodmanson]

Hi Olaf,

First of all, for taking elements of an array to a power in element-by-element fashion, people generally write A.^n instead of power(A,n). The dot indicates element-by-element. If A is a scalar, like you are using here, you don't need the dot: A^n.

I believe you are getting nans because of the line in the ratio function, which I will rewrite as

factor = var^ex / ( K^ex + var^ex) .

ex = n in one of your terms, so ex gets large. Eventually both var^ex and K^ex are driven below the smallest supportable number to zero,

>> realmin  ans = 2.2251e-308

and you get 0/0 = nan. To avoid this, you can use the equivalent

factor = 1 / ( (K/var)^ex + 1 )

The ( .. )^ex term can go to zero or infinity, but in neither case do you get nan.

Incidentally, using 'factor' in the ratio function is all right because the variable is local to the function. But 'factor' is also the name of a Matlab function and it's best not to use that name as a variable. Steering people away from names of existing functions is not one of Matlab's strong points.

[Olaf Sikorski]

Wow, the mistake was simpler than I thought. Many thanks for identifying the error and the comments. I tried it and it solved the issue, so I can continue the actual work.

Thank you very much again!