.. _ode:
Numeric Differential Equations
================================
.. index:: differential equations, ode
.. topic:: Reading Assignment
Please read chapter 8 of *Physical Modeling in MATLAB*, by Allen B.
Downey [DOWNEY11]_.
Basic Idea
-----------
Recall the definition of a derivative:
.. math:: \frac{dy}{dt} = \lim_{h \rightarrow 0} \frac{y(t + h) - y(t)}{h}.
Considered in terms of discrete samples of :math:`y`, the derivative
becomes
.. math:: y_i' = \frac{y_{i+1} - y_i}{h}.
Rearranged to find :math:`y_{n+1}` given the previous value, :math:`y_n`,
and the slope of the function, :math:`y_n'`, we have
.. math:: y_{i+1} = y_i + h\,y_i'
This equation is the basis for the simplest ODE solver, Euler's method.
Differential equations are given in terms of an equation for the derivative,
:math:`y' = f(t, y)` and initial conditions, :math:`y_0`.
.. math:: \spalignsys{y_1 = y_0 + h\,f(t_{0} , y_0);
y_2 = y_1 + h\,f(t_{1} , y_1);
y_3 = y_2 + h\,f(t_{2} , y_2);
\. \vdots;
y_n = y_{n-1} + h\,f(t_{n-1} , y_{n-1})}
Euler's method nicely introduces the concept behind numerical methods to
solve differential equations. However, results from Euler's method lack
accuracy because it assumes a fixed slope over each interval
:math:`\{t_i, t_{i+1}\}`, which is usually not the case.
This `MathWorks web page
`_
has some good tutorial videos about numerical solutions to differential
equations. Please watch at least the first three videos (*Euler, ODE1*,
*Midpoint Method, ODE2*, and *Classical Runge-Kutta, ODE4*). Although, the
ODE solvers in MATLAB are slightly more complex than the examples in the
first three videos, they explain the basic concepts.
As explained in Dr. Moler's videos, the Midpoint Method, Runge-Kutta, and
MATLAB's suite of ODE solvers, all follow the basic concept of Euler's
method, but with more sophisticated methods for calculating the slope at
each step. MATLAB's ODE solvers also adapt the size of :math:`h` for
efficiency.
ode45
-----
The primary ODE solver is called ``ode45``. It has good accuracy and
usually runs fast enough. We set-up to use ``ode45`` by defining a function
handle to define the derivative
.. math:: \frac{dy}{dt} = f(t, y)
Then ``[t,y] = ode45(f,tspan,y0)``, where ``tspan = [t0 tf]``,
integrates the system of differential equations :math:`y'=f(t,y)` from
:math:`t_0` to :math:`t_f` with initial conditions :math:`y_0`. Each row in
the solution array ``y`` corresponds to a value returned in column vector
``t``. Note that the function handle ``f`` should take two arguments,
although in many cases, the function only operates on one variable --
usually the second, ``y``, variable.
When no output is saved, ``ode45`` will plot the function :math:`y(t)`.
The Logistic Equation
----------------------
This example is a model of growth slowed down by competition. The
competition is within the population for limited resources. The growth gives
a :math:`y` term in the ODE -- the growth is proportional to the
population size. The competition (:math:`y` against :math:`y`) adds a
:math:`-y^2` term.
.. math:: \frac{dy}{dt} = a\,y - b\,y^2.
The solution is modeled with positive and negative values for :math:`t`.
The analytic solution is: ([STRANG14]_ page 55)
.. math:: y(t) = \frac{a}{b + d\,e^{-t}}.
Where, :math:`d = \frac{a}{y(0)} - b`.
In the distant past, (:math:`t = -\infty`), the population size is
modeled as zero. The steady state (:math:`t = \infty`) population size
is :math:`y(t) = K = \frac{a}{b}`, which is the sustainable population size.
The half way point, can be modeled as :math:`y(0) = a/2b`.
We will compute the positive half of the solution and then reflect the
values to get the negative time values. I'm not yet aware of how to specify
an initial condition that does not start at the lower value of the ``tspan``
range. Please let me know if you find out how this is done.
::
% Logistic Equation: a = 1, b = 1
f = @(t, y) y - y.^2;
tspan = [0 5];
y0 = 1/2;
[t,y] = ode45(f,tspan,y0);
figure, plot(t,y)
t1 = flip(-t);
y1 = flip(1-y);
t2 = [t1;t];
y2 = [y1;y];
analytic = 1./(1 + exp(-t2));
figure, plot(t2,y2,t2,analytic)
title('Logistic Equation')
legend('ode45 Numeric Solution','Analytic Solution', ...
'Location', 'northwest');
When you run this code, you will see that the numeric solution from
``ode45`` is pretty much identical to the analytic solution.
Relation to Integration
------------------------
Since solving both an integration problem and a differential equation are
reverse operations of taking the derivative of a function, then when the
differential equation does not contain :math:`y` terms there is a
relationship between the value from a definite integral and the difference
between the first and last solution points of the differential equation.
If :math:`\frac{dy}{dx} = f(x)`, then
.. math:: \int_a^b f(x) dx = y(b) - y(a).
For example, consider
.. math:: \int_0^1 \frac{4}{1 + x^2} dx = \pi
Solve by symbolic integration::
>> syms x
>> int(4/(1+x^2))
ans =
4*atan(x)
>> int(4/(1+x^2),0,1)
ans =
pi
Solve by numeric integration::
>> f = @(t) 4./(1 + t.^2);
>> integral(f,0,1)
ans =
3.1416
Solve by numeric differential equation::
>> y0 = 0;
>> eqn = @(t,y) 4./(1 + t.^2);
>> [t, y] =ode45(eqn,[0 1],y0);
>> y(end)
ans =
3.1416
Other ODE Solvers
------------------
MATLAB has several ODE solvers, they are reviewed on this `web page
`_.
As noted on the web page, some ODE are considered *stiff*. When solvers such
as ``ode45`` attempt to solve stiff problems, the step size becomes very
small and the solver runs very slow. If this occurs, try one of the stiff
solvers such as ``ode15s``.
.. seealso::
MATLAB documentation on solving `Systems of ODEs
`_,
and `Higher-Order ODEs
`_.