Solving simultaneous time-dependent matrix equations

I am attempting to solve a set of time-dependent equations which involve matrices. I have tried to use both ode45 and ode15 which both work for equations with one-dimensional variables; as follows:

% Defining constants
A_a = 7.7790;
g1_1 = 0.0019;
g2_1 = 0.0021;
ga_1 = 0.1862;
Hv1 = 2.375;
Hv2 = 2.375;
H1 = 3.9;
CE = 0.82;
Astar1 = 0.9282;
Astar2 = 1.1602;

tspan = [0 9];
Q0 = zeros(3,1);
[t,Q] = ode45(@(t,Q) fun4(t,Q,A_a,g1_1,g2_1,Hv1,Hv2,ga_1,H1,CE,Astar1,Astar2), tspan, Q0)
 plot(t,Q(:,1),'-o',t,Q(:,2), '.-')

function dQdt = fun4(t,Q,A_a,g1_1,g2_1,Hv1,Hv2,ga_1,H1,CE,Astar1,Astar2)
dQdt = zeros(3,1);
dQdt(1) = (Astar1^2)*((g1_1*Hv1 + ga_1*(H1 + CE))-((Q(3)^2)/A_a^2));
dQdt(2) = (Astar2^2)*((g2_1*Hv2 + ga_1*(CE))-((Q(3)^2)/A_a^2));
dQdt(3) = Q(1) + Q(2);
end

This solves the equation as expected. However when the constants g1_1, g2_1, and ga_1 each become a 1-by-2 matrix, I tried using a similar function (given below), however, it will not solve the code. Is there a way to use matrices in differential equations?

% Defining constants
A_a = 7.7790;
g1_1 = [0.0019 0.0076]
g2_1 =  [0.0021 0.0083]
ga_1 =  [0.1862 0.3725]
Hv1 = 2.375;
Hv2 = 2.375;
H1 = 3.9;
CE = 0.82;
Astar1 = 0.9282;
Astar2 = 1.1602;

tspan = [0 9];
Q0 = zeros(3,2);
[t,Q] = ode45(@(t,Q) fun4(t,Q,A_a,g1_1,g2_1,Hv1,Hv2,ga_1,H1,CE,Astar1,Astar2), tspan, Q0)%, options)
 plot(t,Q(:,1),'-o',t,Q(:,2), '.-')

function dQdt = fun4(t,Q,A_a,g1_1,g2_1,Hv1,Hv2,ga_1,H1,CE,Astar1,Astar2)
dQdt = zeros(3,2);
dQdt(1) = (Astar1^2)*((g1_1*Hv1 + ga_1*(H1 + CE))-((Q(3)^2)/A_a^2));
dQdt(2) = (Astar2^2)*((g2_1*Hv2 + ga_1*(CE))-((Q(3)^2)/A_a^2));
dQdt(3) = Q(1) + Q(2);
end

Answer 1

[Stephan]

Hi,

try:

% Defining constants
A_a = 7.7790;
g1_1 = [0.0019 0.0076]
g2_1 = [0.0021 0.0083]
ga_1 = [0.1862 0.3725]
Hv1 = 2.375;
Hv2 = 2.375;
H1 = 3.9;
CE = 0.82;
Astar1 = 0.9282;
Astar2 = 1.1602;
tspan = [0 9];
Q0 = zeros(3,2);
[t,Q] = ode45(@(t,Q) fun4(t,Q,A_a,g1_1,g2_1,Hv1,Hv2,ga_1,H1,CE,Astar1,Astar2), tspan, Q0)%, options)
plot(t,Q(:,1),'-o',t,Q(:,2), '.-')

function dQdt = fun4(t,Q,A_a,g1_1,g2_1,Hv1,Hv2,ga_1,H1,CE,Astar1,Astar2)
dQdt = zeros(6,1);
dQdt(1:2) = (Astar1^2)*((g1_1*Hv1 + ga_1*(H1 + CE))-((Q(3)^2)/A_a^2));
dQdt(3:4) = (Astar2^2)*((g2_1*Hv2 + ga_1*(CE))-((Q(3)^2)/A_a^2));
dQdt(5:6) = Q(1) + Q(2);
end

You need to look if the plot has to be changed, since the solution Q now has 6 columns. One way would be:

subplot(3,1,1)
plot(t,Q(:,1),'-o',t,Q(:,2), '.-')
subplot(3,1,2)
plot(t,Q(:,3),'-o',t,Q(:,4), '.-')
subplot(3,1,3)
plot(t,Q(:,5),'-o',t,Q(:,6), '.-')

Also im not sure about the last line in the function - maybe:

dQdt(5:6) = Q(1:2) + Q(3:4);

?

Best regards

Stephan

[JS]

Thank you for this! This works in calculating each value of Q, however it splits Q1 into two answers. The values of g1_1, g2_1, and ga_1 are also time-dependent and were calculated earlier (in reality each is a 41-by-1 matrix rather than a 2-by-1 matrix, though I was trialling the problem with only a 2-by-1 matrix) . I was hoping to obtain time-dependent answers to Q1, Q2, and Q3 which vary in time as with ga_1, g1_1 and g2_1, and will also be of dimension 41-by-1 which can be plotted against time (time is a 41-by-1 matrix).

Perhaps this could be represented better by using tspan = [0:9:9] in the code above (so that time is also a 2-by-1 matrix).

[Stephan]

your function as you wrote it will give as many answers for Q1 as g1_1 has elements, because it follows the rules of linear algebra. You add and multiply scalar values to a 1x2 vector, which always results in a 1x2 vector.

If i understood you correctly, you should combine the previous time dependent calculation of the factors with the calculation of Q in one common function to get what you want.