ode45 with matrix 1st order ode

1 view (last 30 days)
Haya M
Haya M on 14 Nov 2020
Commented: Haya M on 15 Nov 2020
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?

Accepted Answer

Alan Stevens
Alan Stevens on 15 Nov 2020
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 Comments
Alan Stevens
Alan Stevens on 15 Nov 2020
Edited: Alan Stevens on 15 Nov 2020
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 on 15 Nov 2020
yes, I got it! Thank you very much Alan.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!