The matrix Riccati differential equation:
dX/dt = A'X + XA - XBB'X + Q
function dXdt = mRiccati(t, X, A, B, Q) X = reshape(X, size(A)); %Convert from "n^2"-by-1 to "n"-by-"n" dXdt = A.'*X + X*A - X*B*B.'*X + Q; %Determine derivative dXdt = dXdt(:); %Convert from "n"-by-"n" to "n^2"-by-1
Then, you can use the ODE45 function to solve this problem:
[T X] = ode45(@(t,X)mRiccati(t, X, A, B, Q), [0 10], X0)
For example, using the sample data:
A = [1 1; 2 1]; B = [1; 1]; Q = [2 1; 1 1]; X0 = [1; 1; 1; 1];
You can use the following command to solve the system of differential equations:
[T X] = ode45(@(t,X)mRiccati(t, X, A, B, Q), [0 10], X0)
ODE45 returns "X" as a vector at each time step. You may use the following code to reshape each row of "X" to get the matrix and store it in a cell array:
[m n] = size(X); XX = mat2cell(X, ones(m,1), n); fh_reshape = @(x)reshape(x,size(A)); XX = cellfun(fh_reshape,XX,'UniformOutput',false);
The results of this can be verified by the LQR function:
[K,S,E] = lqr(A, B, Q, 1)
where "S" should have results very similar to the last elements in "X" or "XX". The LQR function computes the steady-state value of the system. In this example, we generated the solution for up to "t = 10", which is an adequate approximation of infinity for this problem.
For more information on ODE45 and other such solvers, refer to the function reference page for ODE45 in the MATLAB documentation.