1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving Systems of ODE's

  1. Mar 22, 2008 #1
    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

    [tex]dx_1 = 2x_2[/tex]
    [tex] dx_2 = -2x_1[/tex]
    [tex] dx_3 = 200x_4[/tex]
    [tex]dx_4 = -200x_3[/tex]

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

    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.
     
    Last edited: Mar 22, 2008
  2. jcsd
  3. Mar 22, 2008 #2
  4. Mar 22, 2008 #3
    It's a linear system, so the solution is
    [tex] Y(t) = e_m^{At}x_0[/tex]
    where x_0 is the initial condition, e_m is the matrix exponential, and A is a the matrix
    [tex]
    \begin{bmatrix}
    0 & 2 & 0 & 0 \\
    -2 & 0 & 0 & 0 \\
    0 & 0 & 0 & 200 \\
    0 & 0 & -200 & 0
    \end{bmatrix}
    [/tex]
    A is defined so that dx = A*x is the ode function in your code.

    The matrix exponential is defined as
    [tex]
    e_m^{A} = T e^{\Lambda} T^{-1}
    [/tex]
    where
    [tex]
    A = T \Lambda T^{-1}
    [/tex]
    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:
    Code (Text):

    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
    Code (Text):

    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.
     
  5. Mar 22, 2008 #4
    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?
     
  6. Mar 22, 2008 #5
    From your Matlab code:
    Code (Text):

    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:
    Code (Text):

    [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
    [tex]\dot{x} = ax[/tex]
    then the solution is an exponential
    [tex]x(t) = x_0 e^{at}[/tex]
    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
    [tex]\dot{x} = Ax[/tex]
    where A is a matrix, then factor it via eigendecomposition to obtain
    [tex]\dot{x} = T^{-1}DTx[/tex]
    Then change variables to [tex]z = Tx[/tex] and left-multiply by T to obtain
    [tex]\dot{z} = Dz[/tex]
    which is now 4 uncoupled equations because D is diagonal! So each one has the solution
    [tex]z_i(t) = z_{i0} e^{D_it}[/tex]
    and
    [tex]z(t) = e^{Dt}z_0[/tex]
    Now change back to x via [tex]x = T^{-1}z[/tex] to obtain
    [tex]x(t) = T^{-1}e^{Dt}Tz_0[/tex]
    which is the solution presented originally. I think I switched T and inverse of T, but it doesn't really matter.

    Matlab:
    Code (Text):

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Solving Systems of ODE's
  1. Solving system of ODE? (Replies: 7)

Loading...