28 views (last 30 days)

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;

else

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));

end

end

time(j+1) = time(j) + dt;

end

% 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

plot(time,Vm_num);

hold on

plot(time,Vm_ana);

title('Comparing Analytical and Numerical Membrane Potentials over Time');

xlabel('Time(ms)');

ylabel('Membrane Potential (mV)');

legend('Numerical Vm','Analytical Vm');

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.

Star Strider
on 6 Apr 2020

The loop is not necessary. Use the vectorisation abilities MATLAB provides, as I did to calculate ‘Vm_ana’.

This revised code makes all the vector and matrix sizes work correctly, however there are obviously still problems that you need to sort:

% 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

lambda2 = 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(fix(N),fix(TI));

time = zeros(1,fix(N));

% 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;

else

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));

end

end

time(i) = (i-1) + dt;

end

% Calculate the steady state of the cable equation

lambda = sqrt(lambda2);

x = linspace(-L/2, L/2, numel(time));

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

plot(time,Vm_num);

hold on

plot(time,Vm_ana);

title('Comparing Analytical and Numerical Membrane Potentials over Time');

xlabel('Time(ms)');

ylabel('Membrane Potential (mV)');

legend('Numerical Vm','Analytical Vm');

We’ll work on ‘lambda’ when I understand how you want to define it. Just now, that’s not obvious.

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

plot(x,Vm_num(:,end));

hold on

plot(x,Vm_ana);

instead of

plot(time,Vm_num);

hold on

plot(time,Vm_ana);

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.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/515591-1d-cable-model#comment_821951

⋮## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/515591-1d-cable-model#comment_821951

## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/515591-1d-cable-model#comment_822342

⋮## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/515591-1d-cable-model#comment_822342

Sign in to comment.