What is an ordinary differential equation? It is an equation which has unknown. The unknown is a function of one variable. It has to have derivative(s). What is an initial value problem (IVP)? The standard form is:
y'(t) = f(t,y), y(t0) = y0.
Note: y(t) is the unknown function that we want to solve. Both y(t) and f(t,y) can be vector functions. f(t,y) itself is not a differential equation. If f(t,y) =f(y), the differential equations is called AUTONOMOUS.
Equilibriums (Steady state solutions) For an autonomous system, the solution of the algebraic equation(s) f (y) = 0 is (are) called the equilibrium(s). Ex: y1' = y1 - y12 - y1 y2 The equilibrium are (0,0), and (1,0).
y1' = - y2 + y1 y2Equilibriums can be classified as stable, unstable, saddle point, stable spirals, unstable spirals, and center (corresponding to a limit cycle usually).
The stability of an equilibrium is determined by the eigenvalues of the Jacobian which the derivative of f(y) at the equilibrium.
If the derivative of f(y) has very different magnitude of the eigenvalues, then the ODE system is called a STIFF system. For a stiff system, the solutions may change very rapidly over some transient period.
Change high order equations to a first order system. Ex: y''' - 2 sin(t) y'' + sin(y') + log(y) = 0, y(0) = 0, y'(0)=1, y''(0)=-1.
Let y1 = y, y2 = y', y3 = y'', then we have
y1' = y2 y1 (0) = 0
y2 ' = y3 y2(0) = 1
y3' = 2 sin(t) y3 - sin( y2) + log(y1) y3 (0) = -1
Euler's method: yi+1 = yi + h f( ti, yi), i = 0, 1, 2, Error ~ C h, first order method. Truncation error of the Euler's method: y(t+h) - yi+1 ~ C h2
Implicit Euler's method: yi+1 = yi + h f( ti, yi+1), i = 0, 1, 2, ....
Improved Euler's method: ypredict = yi + h f( ti, yi),
yi+1 = yi + h [ f( ti, yi) + f( ti, ypredict ) ]/2, i = 0, 1, 2, Error ~ C h2, second order method.
Implicit trapezoidal method (no prediction and correction) yi+1 = yi + h [ f( ti, yi) + f( ti, yi+1 ) ] / 2, i = 0, 1, 2, Error ~ C h2,
Taylor's methods. Use the derivatives of f(t,y) at ti to get high order methods. yi+1 = yi + h y'( ti) + h2 y''(ti) /2 + h3 y'''(ti) /6 + h4 y(4)(ti) /24 + ...
Approximate y(n) using f(t,y) and its partial derivatives at ti.
Second order Taylor's method:
yi+1 = yi + h f( ti , yi ) + h2 ( ft (ti, yi ) + fy (ti, yi ) f( ti , yi ) )/2, i = 0, 1, 2, Error ~ C h2
Third order Taylor's method:
yi+1 = yi + h f( ti , yi ) + h2 ( ft (ti, yi ) + fy (ti, yi ) f( ti , yi ) )/2 +
h3 ( ftt + 2 fty f + fyy f2 + fy ft f + fy2 f )/6, i = 0, 1, 2, Error ~ C h3
all the values are evaluated at ( ti , yi ).
Runge-Kutta Methods: Based on the integral equation, or Taylor expansions. Use the function values of f(t,y) between ti and ti+1 only. There are many Runge-Kutta methods with same order of accuracy. A RK(4) Method:
- Select a few points between ti and ti+1
- Choose prediction(s)
- Form a linear combination to get better approximations.
k1 = f( ti , yi )
k2 = f( ti + h/2 , yi + h k1 /2 )
k3 = f( ti + h/2 , yi + h k2 /2 )
k4 = f( ti + h , yi + h k3 )
yi+1 = yi + h (k1 + 2 k2 + 2 k3 + k4 )/6
Newton's colling, y' = c ( ysur - y(t) ) Falling mass. y' = - g + k u2/m Measuring body's temperatures Pendulum and springs: m y'' + c y' + k y = f(t)
Population, A logistic model predator-prey model: y1' = a y1 - b y1 y2 The prey
y2' = - c y2 + d y1 y2 The predator.
Chemical Reactions:
k1
A + B <-----------------> AB
k2k3
2 A + C <-----------------> A2C
k4
Let y1, y2, y3, y4, and y5, represent the concentrations of the substance of A, B, AB, C, A2C, rewrite the second reaction equation as
k3
A + A + C <-----------------> A2C
k4then the differential equation can be written as
y1' = - k1 y1 y2 + k2 y3 - 2 k3 y12y4 + 2 k4 y5
y2' = - k1 y1 y2 + k2 y3
y3' = k1 y1 y2 - k2 y3
y4' = - k3 y12y4 + k4 y5
y5' = k3 y12y4 - k4 y5Usually the differential equation form chemical reactions are still systems.
Usually we put Matlab commands into a file, e.g., eulerr.m, rk4.m, ... etc. In matlab, we can simply type the file name, e.g., eulerr, rk4, the Matlab will execute the matlab commands in the file line by line. We can add a matlab function by creating a .m file, e.g. , fmass.m, fpop.m, it can have inputs and outputs, e.g., function y = fmass(t,y)
y = -32 + .1*abs(y) + .001*y*y;
The file has to be named as fmass.m. In matlab, we can use it as a function:>> f(2.5, 3)
>> z = fmass(t,y);
>> y(k+1) = y(k) + h* f(t(k),y(k));
In Matlab, use ; to compress the output; use , to get the output on the screen.
In matlab, use diary to record everything into a text fiel.>> diary myrecord
>> eulerr
............
>> diary off
Use save to save variables: >> save myrecord t y -ascii % Save t and y into the file myrecord.
An example to solve the spread of disease with incubation included: S' = - a S I
E' = a S I - c E
I' = c E - b I
R' = b IIn order to use Matlab odesuit to solve the problem, we need to set
y1 = S, y2 = E, y3 = I, y4 = R, and the new system is:
y1' = - a y1 y3
y2' = a y1 y3 - c y2
y3' = c y2 - b y3
y4' = b y3We use a matlab file called spread.m to define the right hand side:
function yp = spread(t, y)Then we can write a Matlab code to input the data, call ode solver, and plot the solutions
global a b c
k=length(y);
yp=zeros(k,1);yp(1) = -a*y(1)*y(3);
yp(2) = a*y(1)*y(3)-c*y(2);
yp(3) = c*y(2)-b*y(3);
yp(4) = b*y(3);
clear all; close all
global a b c
a=1; b=5; c=1;
t0=0; tfinal=5; y0=[8.5 0 1 0];[t y]=ode45('spread',[t0,tfinal],y0);
y1 = y(:,1); y2 = y(:,2); y3 = y(:,3); y4 = y(:,4);
plot(t,y1,t,y2,t,y3,t,y4) % The solution plots.
figure(2); plot3(y(1),y(2),y(3)) % One of phase plot.