MATLAB Answers

1D cable model

8 views (last 30 days)
Kate Heinzman
Kate Heinzman on 5 Apr 2020
Commented: darova on 6 Apr 2020
When I try to pre-allocate the vectors before the for-loop it keeps giving me an error by saying the size inputs must be integers and I am confused as to why it's giving me this error. Then if I get rid of the preallocation, it throws an error in the line Vm_num(1,j) = Vm_num(2,j) - E * dx; by saying the index in position 1 is invalid. How do I fix this? Also in calculating the steady state equation of the cable what would be the x value used there based on the given model parameters?
% Script to numerically calculate changes of Vm(x,t) produced by
% uniform-field electric shock in a 1D monodomain uniform passive cable
% model
clear all
% Model parameters
a = 0.001; % cell radius in cm
L = 1; % cable length in cm
Rm = 5000; % specific membrane resistance in ohm*cm^2
Cm = 0.001; % specific membrane capacitance in mF/cm^2
Ri = 500; % intracellular resistivity in ohm*cm
E = 5; % shock field strength in V/cm
T = 50; % shock duration in ms
dx = 0.01; % integration step in cm
dt = 20; % integration step in us
Tm = Rm * Cm; % time constant in ms
lamda2 = a*Rm / 2*Ri;
% Number of nodes in the cable is equal to the cable length divided by the
% spatial integration step
N = L / dx;
% Time iterations is the shock duration divided by the time integration
% step
TI = T / dt;
% Preallocate numerical Vm and time vector
Vm_num = zeros(N,TI);
time = zeros(1,TI);
% i is the index corresponding to the nodes of the cable
for i = 1:N
% j is the time index
for j = 1:TI
% Boundary Condition
if i == 1
Vm_num(1,j) = Vm_num(2,j) - E * dx;
% Boundary Condition
elseif i == N
Vm_num(N,j) = Vm_num(N-1,j) + E * dx;
Vm_num(i,j+1) = Vm_num(i,j) + (dt/Tm)*((lambda2/dx*dx)*(Vm_num(i+1,j)-2*Vm_num(i,j)+Vm_num(i-1,j))-Vm_num(i,j));
time(j+1) = time(j) + dt;
% Calculate the steady state of the cable equation
lambda = sqrt(lambda2);
Vm_ana = (lambda*E)*(sinh(x???/lambda)/cosh(L/2*lambda));
% Plot the final numerical solution Vm(x) and plot the analytical steady
% state solution of the cable equation for Vm(x) in the same figure
hold on
title('Comparing Analytical and Numerical Membrane Potentials over Time');
ylabel('Membrane Potential (mV)');
legend('Numerical Vm','Analytical Vm');
darova on 6 Apr 2020
Use button for code inserting

Sign in to comment.

Answers (1)

Star Strider
Star Strider on 5 Apr 2020
You are calculating ‘N’ and ‘TI’. They may appear to be integers, however they are not. The solution to that is to force them to be integers with the fix function:
% Preallocate numerical Vm and time vector
Vm_num = zeros(fix(N),fix(TI));
time = zeros(1,fix(TI));
There are other options such as ceil, floor and round, however fix is likely most appropriate here.
Running your code after that correction reveals this one:
lamda2 = a*Rm / 2*Ri;
and changing that to:
lambda2 = a*Rm / 2*Ri;
solves it.
Then an error occurs with this line:
Vm_ana = (lambda*E)*(sinh(x???/lambda)/cosh(L/2*lambda));
and I have no idea what to do with it.
I will be happy to help with it (as well any other problems that may appear) if your code has problems after you correct the ‘Vm_ana’ calculation.
Tommy on 6 Apr 2020
A few notes which I hope are helpful:
1. Don't forget about the units. Keeping the fix() is a good idea, but with a total duration of 50 milliseconds and a step duration of 20 microseconds, you should have 2500 time steps, not 2. See my comment above.
2. "Plot the final numerical solution Vm(x) and plot the analytical steady state solution of the cable equation for Vm(x) in the same figure" makes me think you are supposed to plot Vm_ana over x, as well as the final Vm_num over x. Meaning
hold on
instead of
hold on
3. Currently to fill Vm_num, you pick a given node (i) and loop through all times (j) for that node, then move to the next node. But the value of Vm_num(i,j+1) depends on Vm_num(i+1,j), meaning you need to know the next node's Vm_num (from the previous time). I think it makes more sense to loop through all the nodes at a given time, then move to the next time - in other words, switch the order of the loops.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!