Filling array with for loop, non-singleton dimensions and subscripts

I understand the error message but I don't know what to do :(

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in Trabalho2eii (line 180) v(1,k) = abs((dx*(x2-x1)/(dt*(t2-t1))))*0.1; %[m/s]

Code:
%%Bioelectricidade e Bioelectrofisiologia, MIEBB
% TRABALHO PRÁTICO, 2a parte 
% Katarzyna Więciorek, fc45967; Cristina Mendes, fc48484
% Objetivo: Implementacao em Matlab as equaçoes do cabo necessárias para 
% simular a propagaçao dum potencial de açao, seguindo a descriçao feito 
% no Cap. 6 do Plonsey.

%%Input dado pelo usuário
clear all
close all

prompt = {'Insira o valor da corrente de estímulo Is (em unidades de mA/cm):','Insira o valor do tempo total T (em unidades de ms):','Insira o valor do tempo do início da estimulaçao (em unidades de ms):','Insira o valor da duracao do estimulo (em unidades de ms):','Insira o tamanho do passo temporal (em unidades de ms):'};
dlg_title = 'Estimulaçao';
num_lines = 1;
defaultans = {'2.0','10','2','0.1','0.001'};
answer = inputdlg(prompt,dlg_title,num_lines,defaultans);
Ie = str2double(answer{1}); % um corrente de estímulo Is [uA/cm^2]
T = str2double(answer{2}); % tempo total [ms]
ti = str2double(answer{3}); % inicio [ms]
td = str2double(answer{4}); % duracao [ms]
dt = str2double(answer{5}); %[ms]

% O ciclo para testar os limites das variaveis da entrada
while(1)
    if (Ie < 0) || ((ti+td)>T)
       fprintf(2,'Valores incorrectos! Por favor insira novos valores.\n');
       answer = inputdlg(prompt,dlg_title,num_lines,defaultans);
       Ie = str2double(answer{1}); % um corrente de estímulo Is [uA/cm^2]
       T = str2double(answer{2}); % tempo total [ms]
       ti = str2double(answer{3}); % inicio [ms]
       td = str2double(answer{4}); % duracao [ms]
       dt = str2double(answer{5}); %[ms]

      else
          break;
      end
  end
  %%Inicializacao

% tempo
% y = linspace(x1,x2,n) generates n points. The spacing between the points is dt = (T-0)/(np-1).
np = (T/dt)+1.0; % numero de pontos
t = linspace(0,T,np); %[ms]
tf = ti + td; %fim da estimulacao [ms]

% espaco
lf = 30; % comprimento da fibra [cm]
dx = 0.05; % passo espacial [cm]
nn = (lf/dx)+1.0; %número de nós
x = linspace(0,lf,nn); %[cm]

% constantes
Vr = -60.0; %[mV]
Cm = 1.0; %[uF/cm^2]
gKmax(1:nn,1:np) = 36.0; %[mS/cm^2]
gNamax(1:nn,1:np) = 120.0; %[mS/cm^2]
gL(1:nn,1:np) = 0.3;  %[mS/cm^2]
Ri = 30; %[ohm.cm]
Re = 20; %[ohm.cm]

% Potenciais de Nernst
TC = 6.3; %[Celsius]
%TC = 12.6;
%TC = 18.9;
%TC = 25.0;
[ EK, ENa, EL ] = nernst(TC); %[mV]

% Potenciais [mV]
Vm(1:nn,1:np) = -60.0;
dVm = zeros(nn,np);
vm(1:nn,1:np) = Vm - Vr;

% Alfas e betas [ms^-1]
alfan(1:nn,1:np) = an(vm);
betan(1:nn,1:np) = bn(vm);
alfam(1:nn,1:np) = am(vm);
betam(1:nn,1:np) = bm(vm);
alfah(1:nn,1:np) = ah(vm);
betah(1:nn,1:np) = bh(vm);

% probabilidades
n(1:nn,1:np) = 0.31768;
m(1:nn,1:np) = 0.05293;
h(1:nn,1:np) = 0.59612;
dn = zeros(nn,np);
dm = zeros(nn,np);
dh = zeros(nn,np);

%condutancias [mS/cm^2]
gNa(1:nn,1:np) = gNamax.*(m.^3).*h; %[mS/cm^2]
gK(1:nn,1:np) = gKmax.*(n.^4);

% correntes [uA/cm^2]
Ip=zeros(nn,np); % um corrente de estímulo ip [uA/cm]
INa(1:nn,1:np)= gNa.*(Vm-ENa); 
IK(1:nn,1:np) = gK.*(Vm-EK);
IL(1:nn,1:np)= gL.*(Vm-EL);
Ii = INa + IK + IL;
Im = zeros(nn,np);
Iain = zeros(nn,np); % corrente axial intracellular

% alinea e) ii)
v = zeros(1,20);
ath = zeros(1,10);
k = 1;
j = 1;

%%Calculos

for a = 0.0003:0.00006:0.03; %raio [cm]
cm = Cm*2*pi*a; %Capacidade axial da membrana [uF/cm]
ri = Ri/(pi*a^2); %[ohm/cm]
re = Re/(pi*(3*a^2)); %[ohm/cm]

% mesh ratio
mr = dt./(ri.*cm*1e-3*(dx^2)); %<1fprintf('O mesh ratio para a= %0.3g é %0.3g\n',a,mr);

for Ie = 1.0:0.1:2.0

Ip(1,ti/dt+1:((ti+td)/dt)+1)=-Ie; %primeiro no [mA/cm]
Ip(nn,ti/dt+1:((ti+td)/dt)+1)=Ie; %ultimo no

for i = 1:np-1 % tempo

      %corrente transmembranar (uA/cm^2)
      Im(1,i) = 1000*(1/((pi*a)*(ri + re))).*(((Vm(2,i)-Vm(1,i))/(dx^2))-(re.*Ip(1,i)));
      Im(nn,i) = 1000*(1/((pi*a)*(ri + re))).*(((Vm(nn-1,i)-Vm(nn,i))/(dx^2))-(re.*Ip(nn,i)));    
      Im(2:nn-1,i) = 1000*(1/((2*pi*a)*(ri + re))).*(((Vm(1:nn-2,i)-2*Vm(2:nn-1,i)+Vm(3:nn,i))/(dx^2))-(re.*Ip(2:nn-1,i)));

      %potencial de membrana [mV]
      dVm(:,i)= (dt/Cm)*(Im(:,i)-Ii(:,i));
      Vm(:,i+1) = Vm(:,i) + dVm(:,i);
      vm(:,i+1) = Vm(:,i) - Vr;

      %alfas e betas [ms^-1]
      alfan(:,i) = an(vm(:,i));
      betan(:,i) = bn(vm(:,i));
      alfam(:,i) = am(vm(:,i));
      betam(:,i) = bm(vm(:,i));
      alfah(:,i) = ah(vm(:,i));
      betah(:,i) = bh(vm(:,i));

      %n, m, h
      dn(:,i) = dt.*(alfan(:,i).*(1-n(:,i)) - betan(:,i).*n(:,i));
      dm(:,i) = dt.*(alfam(:,i).*(1-m(:,i)) - betam(:,i).*m(:,i));
      dh(:,i) = dt.*(alfah(:,i).*(1-h(:,i)) - betah(:,i).*h(:,i));
      n(:,i+1)= n(:,i)+dn(:,i);
      m(:,i+1)= m(:,i)+dm(:,i);
      h(:,i+1)= h(:,i)+dh(:,i);

      %condutancias [mS/cm^2]
      gNa(:,i+1) = gNamax(:,i).*((m(:,i+1)).^3).*h(:,i+1); %[mS/cm^2]
      gK(:,i+1) = gKmax(:,i).*((n(:,i+1)).^4);

      % correntes [uA/cm^2]
      INa(:,i+1)= gNa(:,i).*(Vm(:,i)-ENa); 
      IK(:,i+1) = gK(:,i).*(Vm(:,i)-EK);
      IL(:,i+1) = gL(:,i).*(Vm(:,i)-EL);
      Ii(:,i+1) = IK(:,i) + INa(:,i) + IL(:,i);
      %corrente axial intracellular, [mA/cm]
      Iain(1:nn-1,i)=(-1/(ri+re)).*((Vm(2:nn,i)-Vm(1:nn-1,i))/dx-re.*Ip(1:nn-1,i)); 
      Iain(nn,i)=(-1/(ri+re)).*(-Vm(nn,1))/dx-re.*Ip(nn,i); 
  end

if max(max(Vm))>0
ath(1,j)=Ie;
j=j+1;
break
end

end
%Velocidade
t1 = 9000;
t2 = 10000;
x1 = find(Vm(:,t1)==max(Vm(:,t1)));
x2 = find(Vm(:,t2)==max(Vm(:,t2)));
%fprintf('A velocidade da propagacao é igual %0.3g (m/s)\n',v);
v(1,k) = abs((dx*(x2-x1)/(dt*(t2-t1))))*0.1; %[m/s]
k = k+1;
end

%%Gráficos
figure
%Influência do raio da fibra na velocidade de propagação do Potencial
subplot(1,2,1)
plot(a,v(1,:));
title ('Velocidade de Propagação do Potencial de Acção em função do Raio da Fibra')
xlabel ('Raio (cm)')
ylabel ('Velocidade(m/s)')

% Influência do raio da fibra no limiar de estimulacao
subplot(1,2,2)
plot(a,ath(1,:));
title ('Limiar de estimulacao em função do Raio da Fibra')
xlabel ('Raio (cm)')
ylabel ('Intensidade (uA/cm^2)')

Thank you in advance.

Answer 1

[Guillaume]

x1 = find(Vm(:,t1)==max(Vm(:,t1)));

That is odd that x1 ever gets empty. Possibly your matrix contains some NaN and you're using a an older version of matlab that maybe didn't ignore NaNs in max. (NaN is never equal to NaN, so find would never return anything).

In any case, your follow up code assumes that x1 is scalar, which the above line absolutely does not guarantee. If there are several values equal to the max, x1 is going to be a vector. Since you're only interested in one position of the max, simply get the second return value of max which is the index where it is found. So:

[~, x1] = max(Vm(:, t1));
[~, x2] = max(Vm(:, t2));

This is guaranteed to never be empty even if the max is NaN, unless Vm is also empty.

[Katarzyna Wieciorek]

My matlab is 2015b so it is not old. If I have any NaN it is really bad but I don't know how to find them with dbstop because I use functions to fill arrays so it claims NaN is in something like fullpath function (matlab function).

I corrected code as you wrote, it doesn't give any error but velocity is all the time 0 what means it is not calculated. For loop for variable a is probably correct because I am given correct mesh ratio, but two last paragraphs of the code (before plots) are not giving me desired results.

[Guillaume]

I've not really tried to understand your code, I've just ran it with the default values. The values in your Vm matrix start diverging at column 2003, row 601 and then fills up with NaN.

I'm afraid it's for you to find out why. There's obviously an error in one of your formulas. Possibly, it's because of a typo on line 134 where vm is supposed to be Vm (matlab is case sensitive).

If you did really mean to create two variables that differ only by capitalisation (same for cm and Cm on line 112), I would advise you not to do that, it leads too easily to bugs.

Answer 2

[Star Strider]

I did not run your code, however replacing the ‘1’ with a colon ‘:’ in the ‘v’ assignment and vectorising could solve the problem:

v(:,k) = abs((dx*(x2-x1)./(dt*(t2-t1))))*0.1; %[m/s]

[Katarzyna Wieciorek]

Unfortunately:

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in Trabalho2eii (line 179) v(:,k) = abs((dx*(x2-x1)/(dt*(t2-t1))))*0.1; %[m/s]

>> Trabalho2eii Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in Trabalho2eii (line 179) v(:,k) = abs((dx*(x2-x1)./(dt*(t2-t1))))*0.1; %[m/s]

[Star Strider]

I inserted this line just before the one that is throwing the error to see what the result of the calculation is:

QQ = abs((dx*(x2-x1)./(dt*(t2-t1))))*0.1

and the result was:

QQ =
   Empty matrix: 0-by-1

with k=1, because ‘x1’ and ‘x2’ are also empty.

That is the problem. You need to find the reason the calculations yield an empty variable.

[Katarzyna Wieciorek]

x1 = find(Vm(:,t1)==max(Vm(:,t1)));
x2 = find(Vm(:,t2)==max(Vm(:,t2)));

What I am trying to do is to find position, the index of maximum value of Vm in particular column. Is there better way to find it easily?

[Star Strider]

The find function will return all values that are equal to the maximum (there can be more than one), and that is what I thought you wanted so I did not change it.

If you only want the first value and its index, use max with both outputs. The second output is the index.

I would examine the calculations that create the matrix that you then use to create ‘x1’ and ‘x2’ to determine the reason the maximum values return an empty matrix. Look especially for NaN values that are the result of a 0/0 or Inf/Inf operation somewhere in your calculations.