ode45 with matrix 1st order ode

1 visualizzazione (ultimi 30 giorni)
Haya M
Haya M il 14 Nov 2020
Commentato: Haya M il 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?

Risposta accettata

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

Accedi per commentare.

Più risposte (0)

Tag

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by