Infinite Loop, Values not Updating

1 view (last 30 days)
Irene Zhou
Irene Zhou on 26 Nov 2020
Commented: Irene Zhou on 27 Nov 2020
In the while loop, I am comparing the accuracy of Euler's and Midpoint methods (the methods themselves should be correct). If they differ by a certain amount, I reduce the step size by half and redo the calculation for that time point. Else, I update the indices and move on. Plot(t, x(:,4), t, y(:,4)) should gnerate a graph that portrays an action potential.
Somehow after dif hits the specified threshold, x, y, and dif stop to update. The step size keeps getting reduced and I get stuck in an infinite loop. Does anyone have any insights on what might have happened there? Any thought is appreciated!
dt = 0.08;
Vm = 0;
t(1)=0;
i_inj = 17;
%initial values for n, h, m, V
x = [0.317676914060697 0.596120753508460 0.0529324852572496 0]; %set first value to steady state to 0mV
y = [0.317676914060697 0.596120753508460 0.0529324852572496 0];
i = 1;
k = 1;
while t(k)<11
%Euler's method
x(i+1,:) = x(i,:)+dt.*derivs(t(k),x(i,:),i_inj);
%midpoint method
k1 = dt.*derivs(t(k),y(i,:),i_inj);
k2 = dt.*derivs(t(k)+0.5*dt,y(i,:)+0.5.*k1,i_inj);
y(i+1,:) = k2+y(i,:);
dif = abs(x(i+1,4)-y(i+1,4));
if (dif>=0.1)
dt = dt*0.5;
else
t(k+1) = dt+t(k);
i = i+1;
k = k+1;
end
end
function y = derivs(t,x,I_inj)
g_Kbar = 36;
E_K = -12;
g_Nabar = 120;
E_Na = 115;
Gm = 0.3;
V_rest = 10.613;
V=x(4);
if(V == 10)
alpha_n = 1/10;
else
alpha_n = (10-V)./(100*(exp(1-V/10)-1));
end
beta_n = 0.125*exp(-V/80);
n=x(1);
dndt = alpha_n*(1-n)-beta_n*n;
y(1) = dndt;
alpha_h = 0.07*exp(-V/20);
beta_h = 1/(exp(3-V/10)+1);
h=x(2);
dhdt = alpha_h*(1-h)-beta_h*h;
y(2) = dhdt;
if (V==25)
alpha_m = 1;
else
alpha_m = (25-V)/(10*(exp(2.5-V/10)-1));
end
beta_m = 4*exp(-V/18);
m=x(3);
dmdt = alpha_m*(1-m)-beta_m*m;
y(3) = dmdt;
if t>=1 && t<=1.5
dVdt = I_inj+g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
else
dVdt = g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
end
y(4) = dVdt;
end
  2 Comments
KSSV
KSSV on 26 Nov 2020
What is this function derivs?
Irene Zhou
Irene Zhou on 26 Nov 2020
It's a derivative function that calculates and spits out four variables in a row. This function should be correct as I have tested it in other places.
function y = derivs(t,x,I_inj)
g_Kbar = 36;
E_K = -12;
g_Nabar = 120;
E_Na = 115;
Gm = 0.3;
V_rest = 10.613;
V=x(4);
if(V == 10)
alpha_n = 1/10;
else
alpha_n = (10-V)./(100*(exp(1-V/10)-1));
end
beta_n = 0.125*exp(-V/80);
n=x(1);
dndt = alpha_n*(1-n)-beta_n*n;
y(1) = dndt;
alpha_h = 0.07*exp(-V/20);
beta_h = 1/(exp(3-V/10)+1);
h=x(2);
dhdt = alpha_h*(1-h)-beta_h*h;
y(2) = dhdt;
if (V==25)
alpha_m = 1;
else
alpha_m = (25-V)/(10*(exp(2.5-V/10)-1));
end
beta_m = 4*exp(-V/18);
m=x(3);
dmdt = alpha_m*(1-m)-beta_m*m;
y(3) = dmdt;
if t>=1 && t<=1.5
dVdt = I_inj+g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
else
dVdt = g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
end
y(4) = dVdt;
end

Sign in to comment.

Answers (1)

Andy
Andy on 26 Nov 2020
I don't have real values for n0 etc so I put some random numbers in and the resultant t(k) was generally small, much less than the threshold of 11 that you have set for the while loop. I was just wondering if this is the correct variable for the comparison (should you be using k itself). Can you provide valid examples for your starting values?
  4 Comments
Andy
Andy on 27 Nov 2020
It starts going "wrong" on the 33 run through the loop when t=2.3414, I haven't plotted it but what is the function doing at this point?
Irene Zhou
Irene Zhou on 27 Nov 2020
At the end of my Else statement, I update k = k+1 so t(k) is able to reach 11 in theory.
The purpose of comparing the Euler's method to the midpoint method is to take dynamic step sizes. If the difference is too big, take a smaller step. If the difference is small enough, try a bigger step, which I haven't got to.
It starts to go wrong when dif hits the threshold I set, 0.1 in this case. dt keeps getting reduced by half while nothing else gets updated, which is what I am trying to figure out.

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!