solve matrix differential equations with ode45

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

[Torsten]

And what does a matrix as result mean ?

[mpz]

Sorry I forgot to say to use initial condition is the identity matrix.

Essentially, I need to solve the state transition matrix analytically as already done above and numerically as I am trying to do now. So essentially, plotting y for the analytical for tspan = [0 2] will be the same solving after solving it numerically.

[Torsten]

Sorry I forgot to say to use initial condition is the identity matrix.

?

You get a 2x1 vector of functions as solution of your system of ODEs, not a matrix.

[mpz]

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]

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)

Answer 1

[Torsten]

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 = 

[mpz]

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]

Ok. I ran your code and both solutions (at least at first glance) coincide.

So everything is fine now ?