ode45 with matrix 1st order ode

1 回表示 (過去 30 日間)
Haya M
Haya M 2020 年 11 月 14 日
コメント済み: Haya M 2020 年 11 月 15 日
Hi everyone, I'm trying to solve a first order ode in a matrix form using ode45:
where and
$Q=\begin{pmatrix} \sin(x) & 0 \\
0 & \cos(x) \\
\end{pmatrix}$ on .
Here is my code:
clear all
z = 0, %parameter
n = 2;
T0 = eye(n,n)
xspan = [0 5*pi];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[Tl] = ode45(@(x,T) odeTL(x,T,z,n),xspan,T0,opts);
[q,~] = qr(Tl);
Tl = q;
T0 = Tl;
[Tl] = ode45(@(x,T) odeTL(x,T,z,n),xspan,T0,opts);
x1 = 5*pi;
T = deval(Tl,x1);
function dTdx = odeTL(x,T,z,n)
Q = [sin(x) 0;0 cos(x)];
V = Q-z*eye(n,n);
W = T+eye(n,n);
R = T-eye(n,n);
Omg = eye(n,n)-0.5*R'*V*W;
dTdx = T*Omg;
end
As I run the code it said 'Matrix dimensions must agree' and I dont really see how the ode45 works for the matrix case?

採用された回答

Alan Stevens
Alan Stevens 2020 年 11 月 15 日
Does this help
z = -1; %parameter Needs to be negative or T1 blows up.
n = 2;
T0 = eye(n,n);
xspan = [0 5*pi];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[x,Tl] = ode15s(@(x,T) odeTL(x,T,z,n),xspan,T0,opts);
plot(x,Tl),grid
function dTdx = odeTL(x,T,z,n)
T = reshape(T,2,2); % T comes in as a column vector
I = eye(n,n);
Q = [sin(x) 0;0 cos(x)];
V = Q-z*I;
W = T+I;
R = T-I;
O = I-0.5*R'*V*W;
dTdx = T*O;
dTdx = dTdx(:); % dTdx must return as a column vector
end
  3 件のコメント
Alan Stevens
Alan Stevens 2020 年 11 月 15 日
編集済み: Alan Stevens 2020 年 11 月 15 日
Include these lines just before (or just after) the plot command
Tlat5pi = reshape(Tl(end,:),n,n);
disp(Tlat5pi)
The values, of course, will depend on the value you use for z. I chose z = -1 arbitrarily.
Incidentally, I should really have put
T = reshape(T,n,n);
in the loop. However, if you change n to anything other than 2, you will need to specify an equivalent size matrix for Q.
Haya M
Haya M 2020 年 11 月 15 日
yes, I got it! Thank you very much Alan.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

製品

Community Treasure Hunt

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

Start Hunting!

Translated by