How do I solve a system of ODEs with initial conditions in MATLAB?

[ColdFusion85]

For a project we need to do some simulations in matlab. The preliminary part requires us to solve a system of differential equations with initial values given. The system is

dx_1 = 2x_2
dx_2 = -2x_1
dx_3 = 200x_4
dx_4 = -200x_3

x(0) = [1, 0, 1, 0]

My MATLAB code is this:

function part1

options = odeset('RelTol',1e-3,'AbsTol',[1e-3 1e-3 1e-3 1e-3]);
[T,Y] = ode45(@prelim,[0:1:10],[1 0 1 0],options);
[T,Y]

function dx = prelim(t,x)
dx = zeros(4,1); % a column vector
dx(1) = 2*x(2);
dx(2) = -2*x(1);
dx(3) = 200*x(4);
dx(4) = -200*x(3);

end

end

My result is this:

ans =

ans =

0 1.0000 0 1.0000 0
1.0000 -0.4161 -0.9093 0.5098 0.8533
2.0000 -0.6536 0.7568 -0.4680 0.8700
3.0000 0.9602 0.2794 -0.9808 0.0445
4.0000 -0.1455 -0.9894 -0.5387 -0.8140
5.0000 -0.8391 0.5440 0.4193 -0.8749
6.0000 0.8439 0.5366 0.9601 -0.0893
7.0000 0.1367 -0.9906 0.5664 0.7728
8.0000 -0.9577 0.2879 -0.3695 0.8780
9.0000 0.6603 0.7510 -0.9373 0.1339
10.0000 0.4081 -0.9129 -0.5932 -0.7303

To see if this is correct I need to solve the system analytically as well. I've searched the web, but I can't seem to find any good info on how to solve a system of ODE's with initial conditions. I've taken linear algebra, and I'm guessing I have to find the eigenvalues of matrix of coefficients, but I'm not sure what I'd do from there. Could anyone step me through the process, or refer me to a good online resource? Thanks in advance.

 

[Mindscrape]

Yep, you are on the right track.

http://tutorial.math.lamar.edu/Classes/DE/RealEigenvalues.aspx

You probably don't need to use the odeset function unless you are trying to save time, but if it is working then go ahead and keep it.

 

[LeBrad]

It's a linear system, so the solution is
Y(t) = e_m^{At}x_0
where x_0 is the initial condition, e_m is the matrix exponential, and A is a the matrix
<br /> \begin{bmatrix}<br /> 0 &amp; 2 &amp; 0 &amp; 0 \\ <br /> -2 &amp; 0 &amp; 0 &amp; 0 \\ <br /> 0 &amp; 0 &amp; 0 &amp; 200 \\ <br /> 0 &amp; 0 &amp; -200 &amp; 0 <br /> \end{bmatrix}<br />
A is defined so that dx = A*x is the ode function in your code.

The matrix exponential is defined as
<br /> e_m^{A} = T e^{\Lambda} T^{-1}<br />
where
<br /> A = T \Lambda T^{-1}<br />
is the eigenvalue decomposition of A. You can derive this yourself if you do a Taylor series expansion of exp(A) and write A as it's eigenvalue decomposition, and you'll find all the inside T and inv(T) terms cancel, and you get just the expansion of exp on the diagonal matrix, which is just elementwise exp, and then you'll have the two T terms left on the ends.

To implement this in matlab:

time = [0:1:10];
x0 = [0;1;1;0];
A = zeros(4);  A(1,2) = 2; A(2,1) = -2;  A(3,4) = 200; A(4,3) = -200;
for it = 2:length(time)
    Y2(:,it) = expm(A*time(it))*x0;
end

uses the matrix exponential expm, or to do it by hand

time = [0:1:10];
x0 = [0;1;1;0];
A = zeros(4);  A(1,2) = 2; A(2,1) = -2;  A(3,4) = 200; A(4,3) = -200;
[T D] = eig(A);
for it = 2:length(time)
    Y3(:,it) = T*exp(D*time(it))*inv(T)*x0;
end

Note there is a difference between exp(.) and expm(.)

You will get slightly different answers because your error tolerance is loose in your code. If you tighten the tolerance you will get the same answers.

 

[ColdFusion85]

How do I narrow my error tolerance? I have no experience in MATLAB. Also, could you explain what you just wrote above? Is this how I would solve for the eigenvalues, or is this solving the ODE?

 

[LeBrad]

From your Matlab code:

options = odeset('RelTol',1e-3,'AbsTol',[1e-3 1e-3 1e-3 1e-3]);

the relative tolerance and absolute tolerance are set here to 1e-3. To make the tolerance tighter, just change these numbers to something smaller, for example, 1e-8 (for all of them).

Eigenvalue decomposition:

[T D] = eig(A)

This will give you a matrix T filled with the eigenvectors of A, and a diagonal matrix D which contains the eigenvalues of A on the diagonals. These satisfy A*T = T*D

Solving the ODES:
If you have a single linear ode
\dot{x} = ax
then the solution is an exponential
x(t) = x_0 e^{at}
To solve a system of linear ODEs, you need to separate the 4 coupled equations into 4 uncoupled equation, each of which has the exponential solution above. If you have the system
\dot{x} = Ax
where A is a matrix, then factor it via eigendecomposition to obtain
\dot{x} = T^{-1}DTx
Then change variables to z = Tx and left-multiply by T to obtain
\dot{z} = Dz
which is now 4 uncoupled equations because D is diagonal! So each one has the solution
z_i(t) = z_{i0} e^{D_it}
and
z(t) = e^{Dt}z_0
Now change back to x via x = T^{-1}z to obtain
x(t) = T^{-1}e^{Dt}Tz_0
which is the solution presented originally. I think I switched T and inverse of T, but it doesn't really matter.

Matlab:

time = [0:1:10];  % define vector of time values
x0 = [0;1;1;0];  % define vector of initial conditions
A = zeros(4);  A(1,2) = 2; A(2,1) = -2;  % define A matrix
A(3,4) = 200; A(4,3) = -200;
[T D] = eig(A);      % eigenvalue decomposition of A

% solve the system of ODEs by looping over the time points 
%and computing the solution at each time step and storing it 
% as a column in a matrix Y3
for it = 2:length(time)  % loop through the time steps
    Y3(:,it) = T*exp(D*time(it))*inv(T)*x0;  % solve the ODEs at each time point
end

You should be able to paste this block into your code and it will give you the solution stored in Y3.