Top Level

Ask OLC a question

Athena Consulting Homepage

Helpdesk Stock Answers (for Mac/PC questions)

How to solve ORDINARY DIFFERENTIAL EQUATIONS

Matlab can numerically solve Ordinary Differential equations using 2 methods.

        ODE23 uses 2nd and 3rd order Runge-Kutta formulas
and
        ODE45 uses 4th and 5th order Runge-Kutta formulas

What you first need to do is to break your ODE into a system of 1st order
equations. For instance,

        ..    .    2
        X  + AX + K X = 0       <= Damped harmonic oscillator

can be written as,
        .
        Y1 = Y2
        .              2
        Y2 = - A Y2 - K Y1
                 .
(Let Y1=X and Y2=X)

Now, you need to write a matlab function that takes Y1, Y2, and time
                         .      .
as arguments and returns Y1 and Y2.  To do this, create a file using
emacs or ez or your own favorite editor that contains the following:

		 ----------------------------------------

    function [Ydot] = myode(t,Y)

    % Note: Y(1) => Y1 and Y(2) => Y2
    % t is for time.. if your equations do not contain t, you still have
    % to put a "t" in the function declaration above.

    A = 1.3333;    % Whatever values you want to choose
    K= 2.132;

    Ydot(1) = Y(2);
    Ydot(2) = -A*Y(2)-K^2*Y(1);

    Ydot = Ydot';  % This makes Ydot into a column vector

		 ----------------------------------------

Call this file "myode.m" and save it in a place where matlab knows to look
for it (either the directory from which you started matlab or your own
/mit/username/matlab/ directory).

Now, the moment we have been waiting for: let's solve this ODE.  At the
matlab prompt, type:

	 >> [T,Y]=ode45('myode',[0 10],[0 1]);

Note:

        'myode': Name of the function you wrote above
        [0 10] : Integrate from t=0 to t=10 (seconds)
        [0 1]  : Initial conditions for your system.. [1 0]
                 corresponds to system where,

                        X(t=0) = 1
                        .
                        X(t=0) = 0

                 Which you can think of as a pendulum that has
                 been pulled back and is being let go at t=0.
                 (in this case the pendulum is damped by the A term)

To see the results, type:

        >> plot(T,Y)
or, 
        >> plot(T,Y(:,1),'-',T,Y(:,2),'--');

Here, the solid line corresponds to X (displacement) and the dashed line
   .
to X (speed).

Your coefficients to your equations do not have to be constant.  They can
depend on time.  For instance if your oscillator is damped and driven:

        ..    .    2
        X  + AX + K X = B cos(wt)

Your m-file will look something like this:

		 ----------------------------------------

    function [Ydot] = myode2(t,Y)

    % Note: Y(1) => Y1 and Y(2) => Y2
    % t is for time.. if your equations do not contain t, you still have
    % to put a "t" in the function declaration above.

    A = 1.3333;    % Whatever values you want to choose
    K= 2.132;
    B=.2;
    w=2;

    Ydot(1) = Y(2);
    Ydot(2) = -A*Y(2)-K^2*Y(1)+B*cos(w*t);

    Ydot = Ydot';  % This makes Ydot into a column vector

		 ----------------------------------------

You might try experimenting with the coefficients to see what happens.

    >> [T,Y]=ode45('myode2',[0 10], [0 0]); plot(T,Y);