Error using odearguments (line 93); must return a column vector.

10 views (last 30 days)
Maarten Duijn
Maarten Duijn on 17 Dec 2019
Answered: Marcel Jiokeng on 14 May 2022
Hi guys!
I'm trying to run a stimulation of two neural masses in matlab, however the ode is giving me a error code.
I tried to fix it myself but after some trial and error I keep getting the same error.
Down below is both my run script, which runs the ode, and the function itself.
%Firing rate equations for 2 neural masses with only electrical coupling
%& chemical coupling. N=10000
I_1= pi^2;
I_2= -2*pi^2;
params.tau= 10;
params.g= 1;
params.delta= 1;
N= 10000;
V_1= ones(10000,1)*-65;
R_1= zeros(10000,1);
V_2= ones(10000,1)*-65;
R_2= zeros(10000,1);
I= [I_1;I_2];
y0=[V_1,V_2,R_1,R_2];
options = odeset('Events',@QIF_reset)
tstart=0
tfinal= 100
tout=tstart
yout = y0.';
teout = [];
yeout = [];
ieout = [];
for i = 1:5
% Solve until the first terminal event.
[t,y,te,ye,ie] = ode23(@(t,y)QIF(t,y0,I,params,N),[tstart tfinal],y0,options);
if ~ishold
hold on
end
% Accumulate output. This could be passed out as output arguments.
nt = length(t);
tout = [tout; t(2:nt)];
yout = [yout; y(2:nt,:)];
teout = [teout; te]; % Events at tstart are never reported.
yeout = [yeout; ye];
ieout = [ieout; ie];
% Set the new initial conditions, with .9 attenuation.
y0(:,1) = -100;
y0(:,2) = -65;
y0(:,3) = 0;
y0(:,4) = 0;
% A good guess of a valid first timestep is the length of the last valid
% timestep, so use it for faster computation. 'refine' is 4 by default.
tstart = t(nt);
end
plot(tout,yout(:,1))
And here is the function which the run script runs:
function yp= QIF(~,y,I,params,N)
%QIF ordinary differential equation solver for two neurons electrically
%coupled. The input of these neurons are n1= pi^2 and n2= -2*pi^2
I_1= I(1);
I_2=I(2);
v_1= sum(y(1:10000,1))/N;
v_2= sum(y(1:10000,2))/N;
for i=1:N
V_1d(i)= (((y(i,1).^2)+ I_1)+params.g*(v_1-y(i,1)))/params.tau;
V_2d(i)= (((y(i,2).^2)+ I_2)+params.g*(v_2-y(i,2)))/params.tau;;
R_1d(i)= ((params.delta/(params.tau*pi))+ 2*y(i,3)*y(1)-params.g*y(i,3));
R_2d(i)= ((params.delta/(params.tau*pi))+ 2*y(i,4)*y(2)-params.g*y(i,4));
end
yp = [V_1d;V_2d;R_1d;R_2d];

Answers (2)

Steven Lord
Steven Lord on 17 Dec 2019
Well, does your function return a column vector like the error message states it must, or does it return a row vector?
The easiest solution will probably be to reshape yp to be a column vector before returning it from your function.
  2 Comments
Steven Lord
Steven Lord on 17 Dec 2019
yp = [V_1d;V_2d;R_1d;R_2d];
yp = reshape(yp, 4*N, 1);
Be careful about how your initial conditions are ordered, to ensure they're ordered in the same order as this yp vector. Remember that MATLAB will reshape the matrix in column major order.
V_1d = 1:4;
V_2d = 5:8;
R_1d = 9:12;
R_2d = 13:16;
yp = [V_1d;V_2d;R_1d;R_2d];
yp = reshape(yp, 4*N, 1) % This is not 1:16
If you expected yp to be 1:16, you probably want to preallocate each of your vectors to be column vectors before you fill them in.
V_1d = zeros(N, 1); % and the same for V_2d, R_1d, R_2d
If you preallocate them to be column vectors, then you won't need to use reshape.
V_1d = zeros(N, 1);
V_1d(1) = 1;
V_1d(2) = 4;
V_1d(3) = 9;
V_1d(4) = 16;
V_2d = V_1d + 1;
R_1d = V_1d + 2;
R_2d = V_1d + 3;
yp = [V_1d;V_2d;R_1d;R_2d]

Sign in to comment.


Marcel Jiokeng
Marcel Jiokeng on 14 May 2022
You should initialize yp in the function QIF(~,y,I,params,N) as it follows:
function yp = QIF(~,y,I,params,N)
yp = zeros((with the right length),1);
....
....
....

Community Treasure Hunt

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

Start Hunting!