1D Cable Model

2 views (last 30 days)
Kate Heinzman
Kate Heinzman on 12 Apr 2020
Commented: darova on 13 Apr 2020
Taking the code below how do I fix the plote of Vm(x) at 11 time points t = 0,2,4,...20? I tried doing it as t = [0:2:20], but it kept giving me an error. Then how do I plot curves of Vm(t) at 11 equally spaced x points along the cable. The code I have for this is also giving me errors. Then once that is resolved how do I make a plot of time constants of Vm(t) as a time to reach (1 - (1/e)) of the final Vm value at 11 points along the cable? I have how to find time constants of Vm(t), but I'm not sure how to plot multiple at 11 points along the cable.
% Model parameters
a = 0.001; % cable radius, cm
L = 1; % cable length, cm
Rm = 5000; % membrane resistivity, Ohm*cm^2
Cm = 1; % membrane capacitance, uF/cm^2
Ri = 500; % intracellular resistivity, Ohm*cm
E = 5; % shock field, V/cm
T = 50; % duration of integration, ms
dt = 0.020; % time integration step, ms
dx = 0.0100;% space integration step, cm
tau = Rm*Cm*0.001; % membrane time constant, ms
lambda = sqrt((a*Rm)/(2*Ri)); % space constant, cm
nt = (T/dt) + 1; % number of time points
nx = round(L/dx)+1; % number of nodes
x = -L/2 : dx : L/2; % coordinates of nodes
mesh = dt*lambda^2/(tau*dx^2); % stability requirement mesh < 0.5
% Simplifying equation mathematics for loop
kt = dt/tau;
kx = lambda^2/dx^2;
% Preallocation
vm = zeros(nt,nx);
time = zeros(1,nt);
% Start loop at 2 because of initial condition vm(1,i) = 0
for j = 2 : nt % time loop
for i = 2 : (nx-1) % space loop for internal nodes; Euler method
vm(j,i) = vm(j-1,i) + kt*(kx*(vm(j-1,i+1) - 2*vm(j-1,i) + vm(j-1,i-1)) - vm(j-1,i));
end
% Updated boundary conditions
vm(j,1) = (4*vm(j,2) - vm(j,3) - 2* E*dx)/3; % vm at left edge from boundary conditions
vm(j,nx) = (4*vm(j,nx-1) - vm(j,nx-2) + 2*E*dx)/3; % vm at right edge from boundary conditions
time(j+1) = time(j) + dt;
end
vm = 1000*vm; % convert from V to mV
% Analytical solution
vmAna = 1000*E*lambda*(sinh(x/lambda)/cosh(L/(2*lambda)));
% Plot curves Vm(x) at 11 time points where t = 0,2,4,...20 ms
% Gives an error if try to do t as t = [0:2:20];
t = [2:2:20];
plot(x,vm(t,:));
xlabel('x (cm)');
ylabel('Vm (mV)');
title('Figure 2: L = 1 cm');
% title('Figure 2: L = 0.1 cm');
% Plot curves Vm(t) at 11 equally spaced x points along the cable
xx = linspace(-0.5,0.5,11);
plot(time,vm(:,xx));
xlabel('t (ms)');
ylabel('Vm (mV)');
title('Figure 2: L = 1 cm');
% title('Figure 2: L = 0.1 cm');
% Plot time constants of Vm(t) change as time to reach (1-1/e) of the final
% Vm value at 11 points along the cable
% find time constant of Vm(t) as a time to reach (1-1/e)*Vm_final
j = 1;
level = 1-1/exp(1);
while (abs(vm(j)) < level*vm(nt))
j = j+1;
end;
taunum = j*dt;
  5 Comments
Kate Heinzman
Kate Heinzman on 13 Apr 2020
Ok thanks I fixed it. Can you please answer my other two questions?
darova
darova on 13 Apr 2020
If i understood you correctly: you want to find point where t=1-1/e
Solution
imax = find(t>1-1/exp(1),1,'first'); % index
Then to some plots or something
plot(t(imax),vm(imax,:),'or')

Sign in to comment.

Answers (0)

Categories

Find more on Line Plots in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!