ode45 needs to return a column vector

I am trying to solve the Riccati equation below:

Q = [1 0;
    0 1/2];
R = 1/4
A = [0 1;
    2 -1];
B = [0;
    1];
tspan = [5 0];
P0 = 0;
P = ode45(@(t,P) odefun(t,P,Q,A,B,R), tspan, P0);
function dPdt = odefun(t,P,Q,A,B,R)
dPdt = -Q - A.'*P - P*A + P*B*R.^(-1)*B.'*P;
end

But I am getting the error of

Error using odearguments (line 93)

@(T,P)ODEFUN(T,P,Q,A,B,R) must return a column vector.

Error in ode45 (line 106)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in problem_2 (line 11)

P = ode45(@(t,P) odefun(t,P,Q,A,B,R), tspan, P0);

Answer 1

[Walter Roberson]

Q = [1 0;
    0 1/2];
R = 1/4

R = 0.2500

A = [0 1;
    2 -1];
B = [0;
    1];
tspan = [5 0];
P0 = 0;
P = ode45(@(t,P) odefun(t,P,Q,A,B,R), tspan, P0);

ans = 1×2

2 2

Error using odearguments
@(T,P)ODEFUN(T,P,Q,A,B,R) must return a column vector.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

function dPdt = odefun(t,P,Q,A,B,R)
dPdt = -Q - A.'*P - P*A + P*B*R.^(-1)*B.'*P;
size(dPdt)
end

Q is 2 x 2, so addition or subtraction of Q and something else is going to give you at least a 2 x 2

A is 2 x 2, P is a scalar (because P0 is a scalar), so A.'*P is going to be 2 x 2.

P is a scalar, A is 2 x 2, so P*A is going to be 2 x 2

P is a scalar, B is 2 x 1 sp P*B is 2 x 1. R is a scalar, so P*B*R.^(-1) is 2 x 1. B is 2 x 1 so B.' is 1 x 2, and (2 x 1) * (1 x 2) is 2 x 2. P is scalar, so 2 x 2 * scalar is 2 x 2.

Thus you have multiple terms all of which are 2 x 2.

Your ode function must return a column vector, and the column vector must have the same number of elements as your P0 -- so MATLAB is expecting your function to return a scalar, not 2 x 2. Even if you reshape the 2 x 2 into a column vector, that would give you 4 elements, which would not match the single boundary conditon P0 that you have.

[Tong Chen]

If I change my P0 into a [0 0; 0 0] matrix, will that help?

[Walter Roberson]

Q = [1 0;
    0 1/2];
R = 1/4

R = 0.2500

A = [0 1;
    2 -1];
B = [0;
    1];
tspan = [5 0];
P0 = zeros(2,2);
[t, P] = ode45(@(t,P) odefun(t,P,Q,A,B,R), tspan, P0);
plot(t, P)
legend({'1', '2', '3', '4'})

function dPdt = odefun(t,P,Q,A,B,R)
    P = reshape(P, 2, 2);
    dPdt = -Q - A.'*P - P*A + P*B*R.^(-1)*B.'*P;
    dPdt = dPdt(:);
end

[Tong Chen]

Thank you so much for the help!