ODE system.

  Paradisal Condition:     y'(t) = a y(t).

  It may be valid for a short period, the solution approaches infinity as t approaches infinity.
 

  Add death rate:     y'(t) = a y(t) - b y(t) = (a-b) y(t).

   It  will valishes if b>a, and approaches to infinity if a>b; and remains as a constant if a=b.

   Better model:   y'(t) = a y(t) - b y(t) ^2.
 

 y'(t) = a y(t) - b y(t) z(t)

 z'(t) = - c z(t) + d y(t) z(t).
 

y'''(t) - 2 y''(t) + t y'(t) + 2 y = 0,      y(0) = 0,   y'(0)=1,    y''(0)=0.
 

y₁ = y,
y₂ = y'
y₃ = y''.
 

y₁ ' = y2,
y₂' = y'' = y_3,
y₃' = y''' = 2 y'' - t y' - 2 y = 2 y₃ - t y₂ - 2 y₁
y₁ (0) = y(0) = 0;   y₂(0) = y'(0) = 1;   y₃(0) = y''(0) = 0.

---

 ODE system II

   Format for defining an ODE system:

  function  yp = my_fun(x,y)
     k = length(y);
     yp = zeros(k,1);

     yp(1) =
     yp(2) =

What should be the name of this matlab function?
 
How do you solve the original problem using Matlab?

•       Directly typing.
•       Write a script file,  e.g.  lesson9.m.
 
• Data Structure of the output: (my_euler.m,  and Matlab built-in functions: ode23, ode45, ode15s, ...)

     t                                              y1                 y2                       y3          .............................................

t(1)=t0                                   y1(t0)            y2(t0)                y3(t0)      .............................................

t(2) =t0+h                           y1(t(2))          y2(t(2))             y3(t(2))    .............................................

 .................................................................................................................................................................

t(n+1) = tfinal                     y1(t(n+1))     y2(t(n+1))       y3(t(n+1))   ...........................................
 
 

• How to extract solution components?

 
      y1=y(:,1);  y2=y(:,2);   ................
 
• How to plot solutions against independent variable?

 
                   plot(t,y1),  plot(t,y2)   or
                   figure(2);  plot(t,y2)
 
• How to label plots?         xlabel('t');        ylabel('y_1');       title('Title of the plot')

 

  Runge-Kutta Methods:       One step multi-stage method.

Very popular and useful

More accurate,    h², h³  h⁴ .,....

Adaptive time steps,  different  h.

More stable (later)

 

Key idea:                    y' = f(x,y);  y(x₀) = y₀;                It is equivalent  to:

                  y(x) =  y₀+ int   f(t, y(t))   from    x0    to    x .

From one step to the next,  find  better approximation to the integral.

Try different strategies of integration we get different RK methods.

Examples:    Left  rectangle,  Mid-rectangle,  Right rectangle,   Trapezoidal .....  (class exercise)
 

Steps:  (three groups of coefficients)

• Choose a few points,
• Choose predictions
• Choose final combinations

How to determine the coefficients:

•     Points between  x_i  and x_i+1  are given
•     Other coefficients are determined from the Taylor expansion.

 

Generally, RK(k) method is k-th order accurate, which means that the error is propotional to h^k.

 

 High order methods require more function evaluations (f(x,t)). Therefore if the function is easy to compute, high order methods are more efficient. If the function is too complicated and take too much time to evaluate, then the lower order method such as Euler's method may be more efficient.

Format calling ODE45 (RK(4) + RK(5),  and adaptive time step):

          [x, y] = ode45('yfun',[x0, xfinal] ,y0)

Comparison of different method for the following problem:

    y₁' = -y₂,
    y₂ ' =  y₁
    y₁(0) = 1;  y₂(0) = 0.

The exact solution is   y₁= cos (x),    y₂  =  sin (x),     y₁ + y₂ ' =  1, so phase plot   plot(y1,y2) is a circle.