# ode45 with matrix initial conditions

EDIT: We have a coupled system of 10 ode each. The coupling presents in the last equation. I thought about using a matrix 10 by 2 as initial conditions.

I also followed a similar question with the same title here, but I still get the same errors ('Too many input arguments.')

time = [0 5];
x1_0 = [0 0 0 0 0 0 0 0 0 0];
x2_0 = [0 0 0 0 0 0 0 0 0 0];
initial = [x1_0;x2_0];
x = NaN(length(initial),2*length(time));

[t,x] = ode45(@ode,time,initial);

function [dxdt] = ode(x)
N = 2;
dxdt = NaN(10,2);
A = 3.25e-3;
B = 22e-3;
a = 100;
b = 50;
C = 135;
C1 = C;
C2 = 0.8*C;
C3 = 0.25*C;
C4 = 0.25*C;
C12 = C;

for i = 1:N
dxdt(1,i) = x(6,i);
dxdt(6,i) = A*a*C1*sigm((x(3,i)-x(4,i)+x(5,i))) - 2*a*x(6,i) -
x(1,i)*a^2;
dxdt(2,i) = x(7,i);
dxdt(7,i) = A*a*C2*sigm((x(3,i)-x(4,i)+x(5,i))) - 2*a*x(7,i) -
x(2,i)*a^2;
dxdt(3,i) = x(8,i);
dxdt(8,i) = A*a*C3*sigm((x(1,i))) - 2*a*x(8,i) - x(3,i)*a^2;
dxdt(4,i) = x(9,i);
dxdt(9,i) = B*b*C4*sigm(x(2,i)) - 2*b*x(9,i) - x(4,i)*b^2;
dxdt(5,i) = x(10,i);
if i == 1
j = 2;
elseif i == 2
j = 1;
end
dxdt(10,i) = A*a*C12*sigm(x(3,j)-x(4,j)+x(5,j)) - x(10,i) - x(5,i);
end
end

function X = sigm(u)
u0 = 6e-3;
e0 = 2.5;
r = 0.56e3;
X = 2*e0/(1+exp(r*(u0-u)));
end


If my mistake is the use of a matrix initial conditions instead of a vector, using a 1 by 20 vector, and adjust the ode form accordingly would be not practical, I think - What would be another more efficient way to address the initial conditions -What is unnecessary from the inputs I have given and why? -Is there any other computational way to represent the coupling?

You will need to write your problem such that the unknowns are a single vector, not a matrix. In your example with $$N=2$$, you will have an unknown vector $$x(t)$$ of size $$20\times 1$$ (not a matrix of $$10\times 2$$). You will solve a problem of the following shape $$\dot{x}(t) = A(t)\, x(t), \mathrm{ with }\ x \in \mathbb{R}^{n},\ A(t) \in \mathbb{R}^{n\times n}$$ where $$n = 190 \times 10$$ in your particular case.

If your example is really representative of your problem, the matrix $$A(t)$$ will have constant entries so $$A(t) \equiv A(0) \equiv A$$. It can thus be constructed beforehand (outside of the integration loop). In the function defining the right-hand-side (ode(x) in your case), you will only have a matrix-vector multiplication $$A\, x(t)$$. And given your example, your matrix $$A$$ will be sparse, so you could profit from that using sparse matrix-vector multiplication.

PS. Please note that I think that ode45 expects a function with signature (t,x) as the right-hand-side. So it might be that you have to (formally) code your function as

function dxdt = ode(t,x)
dxdt = A*x
end


Working (from the point of view of syntax, not of simulation) code below. I just linearized your matrix $$m\times n$$ with indices $$(i,j)$$ as a vector with indiceds $$(i+(j-1)*n)$$; reordered the functions, used the correct signature for ode45.

time = [0 5];
%x1_0 = [0 0 0 0 0 0 0 0 0 0];
%x2_0 = [0 0 0 0 0 0 0 0 0 0];
N = 2;
initial = zeros(10*N,1);
%x = NaN(length(initial),2*length(time));

function X = sigm(u)
u0 = 6e-3;
e0 = 2.5;
r = 0.56e3;
X = 2*e0/(1+exp(r*(u0-u)));
end

function dxdt = rhs(t,x)
N=2;
dxdt = NaN(10*N,1);
A = 3.25e-3;
B = 22e-3;
a = 100;
b = 50;
C = 135;
C1 = C;
C2 = 0.8*C;
C3 = 0.25*C;
C4 = 0.25*C;
C12 = C;

for i = 1:N
dxdt(1+(i-1)*10) = x(6+(i-1)*10);
dxdt(6+(i-1)*10) = A*a*C1*sigm((x(3+(i-1)*10)-x(4+(i-1)*10)+x(5+(i-1)*10))) - 2*a*x(6+(i-1)*10) - x(1+(i-1)*10)*a^2;
dxdt(2+(i-1)*10) = x(7+(i-1)*10);
dxdt(7+(i-1)*10) = A*a*C2*sigm((x(3+(i-1)*10)-x(4+(i-1)*10)+x(5+(i-1)*10))) - 2*a*x(7+(i-1)*10) - x(2+(i-1)*10)*a^2;
dxdt(3+(i-1)*10) = x(8+(i-1)*10);
dxdt(8+(i-1)*10) = A*a*C3*sigm((x(1+(i-1)*10))) - 2*a*x(8+(i-1)*10) - x(3+(i-1)*10)*a^2;
dxdt(4+(i-1)*10) = x(9+(i-1)*10);
dxdt(9+(i-1)*10) = B*b*C4*sigm(x(2+(i-1)*10)) - 2*b*x(9+(i-1)*10) - x(4+(i-1)*10)*b^2;
dxdt(5+(i-1)*10) = x(10+(i-1)*10);
if i == 1
j = 2;
elseif i == 2
j = 1;
end
dxdt(10+(i-1)*10) = A*a*C12*sigm(x(3+(j-1)*10)-x(4+(j-1)*10)+x(5+(j-1)*10)) - x(10+(i-1)*10) - x(5+(i-1)*10);
end
end

[t,x] = ode45(@rhs,time,initial);

• You are right with what you suggested and the false is mine. In order to make the post a bit easier to follow through the code, I had to subtract some parts of the system. I couldn't do what you suggested although it would be perfect mach for the way I set the question, because in the original system, there are sigmoid functions that render the system highly nonlinear. This is why I oriented myself towards asking for the appropriate way to set initials – Alex Sep 10 at 12:16
• The actual function is as i edited above. – Alex Sep 10 at 13:36

If this doesn't need to be incredibly fast, just use the "reshape" function to change between matrix and vector representations of the state. This will cause more memory allocations, which can slow things down a lot.

• hm, the 2 couples systems is the first step. Eventually I am going to need to couple 190 of them according to a connectivity matrix and thus, the code needs to be as "light" as it gets – Alex Sep 9 at 10:39
• reshape usually only changes the access logic, the view of the data. The matrix data will still reside in the same flat array. If at all, there should only be minimal memory allocations for a new view instance involved. – LutzL Sep 9 at 16:57