solve matrix differential equations with ode45

12 visualizaciones (últimos 30 días)
mpz
mpz el 27 de Nov. de 2022
Editada: Torsten el 28 de Nov. de 2022
Hi
I want to solve dy/dt=Ay(t) in matlab numerically but having a hard time coming up with a working code. My numerican and analytical results should match. Any help appreciated.
% analytical
A = [0 1;8 -2]
A = 2×2
0 1 8 -2
syms t;
[M,J] = jordan(A);
y = M * expm(J*t) * inv(M)
y = 
% numerical
  5 comentarios
mpz
mpz el 28 de Nov. de 2022
The identity matrix should be 2 by 2 i.e. y(0) = I.
The analytical solution gave me a 4 by 4 matrix which I had to plot y(1,1), y(1,2), y(2,1), y(2,2). I thought that was supposed to be the case solving for dy/dt=Ay(t). Can you please explain why its 2x1 instead. Plotting it will give only 2 plots for y(1,1) and y(2,1). See attached plot I had from the analytical solution
Torsten
Torsten el 28 de Nov. de 2022
But you have to clarify in advance that y is supposed to be a matrix (2x2) matrix of functions, not a 2x1 vector of functions which is the usual interpretation when an ODE is written as
dy/dt = A*y(t)

Iniciar sesión para comentar.

Respuestas (1)

Torsten
Torsten el 27 de Nov. de 2022
Editada: Torsten el 28 de Nov. de 2022
syms t y1(t) y2(t)
A = [0 1;8 -2];
Y = dsolve(diff([y1;y2],t)==A*[y1;y2])
Y = struct with fields:
y2: C1*exp(2*t) + C2*exp(-4*t) y1: (C1*exp(2*t))/2 - (C2*exp(-4*t))/4
y1 = Y.y1;
y2 = Y.y2;
diff(y1,t)-y2
ans = 
0
diff(y2,t)-8*y1+2*y2
ans = 
0
To compare with a numerical solution, you need to prescribe initial conditions for y1 and y2.
Here is your corrected solution:
syms t
A = [0 1;8 -2];
[V,J] = jordan(A)
V = 2×2
-0.2500 0.5000 1.0000 1.0000
J = 2×2
-4 0 0 2
y = V*exp(diag(J)*t)
y = 
  2 comentarios
mpz
mpz el 28 de Nov. de 2022
Editada: Torsten el 28 de Nov. de 2022
Here is how I went about it to solve numerically.
%% Call function file
clear all;clc;close all
syms t;
A = [0 1; 8 -2];
X0 = eye(2);
[T X] = ode45(@(t,X)mRiccati(t, X, A), [0 2], X0);
figure(1)
plot(T,X(:,1),T,X(:,2),T,X(:,3),T,X(:,4))
legend
Here is how I went about it to solve analytical. y(t)= M*e^Jt*M^-1
% analytical
A = [0 1;8 -2]
A = 2×2
0 1 8 -2
syms t;
[M,J] = jordan(A);
phi_s=M * expm(J*t) * inv(M)
phi_s = 
figure(2)
fplot([phi_s(1,1) phi_s(1,2) ...
phi_s(2,1) phi_s(2,2)],[0 2])
%% function file
% f
function dXdt = mRiccati(t, X, A)
X = reshape(X, size(A)); %Convert from "n^2"-by-1 to "n"-by-"n"
dXdt = A.'*X; %Determine derivative
dXdt = dXdt(:); %Convert from "n"-by-"n" to "n^2"-by-1
end
Torsten
Torsten el 28 de Nov. de 2022
Editada: Torsten el 28 de Nov. de 2022
Ok. I ran your code and both solutions (at least at first glance) coincide.
So everything is fine now ?

Iniciar sesión para comentar.

Community Treasure Hunt

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

Start Hunting!

Translated by