Function for ODE45 solver doesn't work
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
First function .m file:
function dVm=dVm(t,y)
% Parameters are time t, y for Vm, m, n and h, input for Istim)
% and given Constants
gNa = 120; % mS/cm^2
gK = 36; % mS/cm^2
gL = 0.3; % mS/cm^2
ENa = 115; % mV
EK = -12; % mV
EL = 10.6; % mV
Cm = 1; % uF/cm^2
% System Functions
dVm(1,1)= @(Istim_sys)(-(1/Cm)*(gNa*((y(2)^3)*y(4)*(y(1)-ENa))+(gK*((y(3)^4)*(y(1)-EK))+(gL*(y(1)-EL))-Istim_sys)));
dVm(2,1) = (0.01*(0.5-y(1)) / (1-exp((9.5-y(1))/10)))*(1-y(2)) - (0.125*exp(-y(1)/ 80))*y(2); % dn/dt where y(2) = n
dVm(3,1) = (0.1*(25-y(1)) / (1-exp((25-y(1))/10)))*(1-y(3)) - (4*exp(-y(1)/ 18))*y(3); % dm/dt where y(3) = m
dVm(4,1) = (0.07*exp(-y(1) / 20))*(1-y(4)) - (1 / (1+exp((30-y(1))/10)))*y(4); % dh/dt where y(4) = h
end
2nd function .m file:
function Istim_sys
Ts1 = 0.5;
Ts2 = 0.55;
Tst = 2;
ts = (0.00:0.001:Ts1); % msec
tsts = ((Ts1+0.001):0.001:Ts2);
tststs = ((Ts2+0.001):0.001:Tst);
t_sys = [ts tsts tststs];
Istim_1 = zeros(1,length(ts));
Istim_2 = 200.* ones(1,length(tsts));
Istim_12 = vertcat(Istim_1',Istim_2');
Istim_3 = zeros(1,length(tststs));
Istim_sys = vertcat(Istim_12,Istim_3');
end
My .m that I build, calling function above:
% Hodgkin-Huxley Model solved with ODE45 MATLAB function
tspan = [0:0.001:2];
% Initial Conditions t0 = 0 msec
Vm0 = 0; % mV
alpha_n0 = 0.01*(9.5-Vm0) / (-1+exp((9.5-Vm0)*0.1));
beta_n0 = 0.125*exp(-(Vm0)/80);
nss0 = alpha_n0 / (alpha_n0+beta_n0); % nss(0mV)
alpha_m0 = 0.1*(25-Vm0) / (-1+exp(0.1*(25-Vm0)));
beta_m0 = 4*exp(-(Vm0)/18);
mss0 = alpha_m0 / (alpha_m0+beta_m0); % mss(0mV)
alpha_h0 = 0.07*exp(-(Vm0)/20);
beta_h0 = 1 / (1+exp((30-Vm0)*0.1));
hss0 = alpha_h0 / (alpha_h0+beta_h0); % hss(0mV)
y0 = [Vm0; nss0; mss0; hss0];
% ODE45 plot
[t,y]=ode45(@dVm,tspan,y0);
% HH parameters
Vm = y(:,1);
n = y(:,2);
m = y(:,3);
h = y(:,4);
% a) Plot Voltage (sol of dV/dt) as function of time
plot(t,Vm);
title('Hodgkin-Huxley Model');
xlabel('Time (ms)');
ylabel('Potential Voltage (mV)');
Error:
Conversion to function_handle from double is not possible.
Error in dVm (line 14)
dVm(2,1) = (0.01*(0.5-y(1)) / (1-exp((9.5-y(1))/10)))*(1-y(2)) - (0.125*exp(-y(1)/ 80))*y(2); % dn/dt where y(2) = n
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in HHproject_a_b_c (line 18)
[t,y]=ode45(@dVm,tspan,y0);
0 Commenti
Risposte (1)
Star Strider
il 16 Nov 2019
You have not provided enough information to determine what the problem is. However, considering that ‘Istim_sys’ may be a vector that may be time-dependent, see the documentation section on ODE with Time-Dependnet Terms as a possible solution.
This is purely a guess, however my interpretation of the problem would be compatible with the error your code has thrown.
4 Commenti
Star Strider
il 16 Nov 2019
The input argument would appear to be ‘t’. I have no idea what the output would be, since I do not understand the function you wrote. You may need to use interp1 to calculate the output.
Vedere anche
Categorie
Scopri di più su Ordinary Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!