Solving simultaneous time-dependent matrix equations

8 次查看（过去 30 天）
JS
JS 2019-2-19
编辑： Stephan 2019-2-19
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

请先登录，再回答此问题。

采纳的回答

Stephan
Stephan 2019-2-19
编辑：Stephan 2019-2-19
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
  2 个评论
JS
JS 2019-2-19
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
Stephan 2019-2-19
编辑：Stephan 2019-2-19
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by