7.5. Numeric Differential Equations

Reading Assignment

Please read chapter 8 of Physical Modeling in MATLAB, by Allen B. Downey [DOWNEY11].

7.5.1. Basic Idea

Recall the definition of a derivative:

\frac{dy}{dt} = \lim_{h \rightarrow 0} \frac{y(t + h) - y(t)}{h}.

Considered in terms of discrete samples of y, the derivative becomes

y_i' = \frac{y_{i+1} - y_i}{h}.

Rearranged to find y_{n+1} given the previous value, y_n, and the slope of the function, y_n', we have

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, y' = f(t, y) and initial conditions, y_0.

\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 \{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 h for efficiency.

7.5.2. 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

\frac{dy}{dt} = f(t, y)

Then [t,y] = ode45(f,tspan,y0), where tspan = [t0 tf], integrates the system of differential equations y'=f(t,y) from t_0 to t_f with initial conditions 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 y(t).

7.5.3. 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 y term in the ODE – the growth is proportional to the population size. The competition (y against y) adds a -y^2 term.

\frac{dy}{dt} = a\,y - b\,y^2.

The solution is modeled with positive and negative values for t. The analytic solution is: ([STRANG14] page 55)

y(t) = \frac{a}{b + d\,e^{-t}}.

Where, d = \frac{a}{y(0)} - b.

In the distant past, (t = -\infty), the population size is modeled as zero. The steady state (t = \infty) population size is y(t) = K = \frac{a}{b}, which is the sustainable population size.

The half way point, can be modeled as 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.

7.5.4. 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 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 \frac{dy}{dx} = f(x), then

\int_a^b f(x) dx = y(b) - y(a).

For example, consider

\int_0^1 \frac{4}{1 + x^2} dx = \pi

Solve by symbolic integration:

>> syms x
>> int(4/(1+x^2))
ans =

>> int(4/(1+x^2),0,1)
ans =

Solve by numeric integration:

>> f = @(t) 4./(1 + t.^2);
>> integral(f,0,1)
ans =

Solve by numeric differential equation:

>> y0 = 0;
>> eqn = @(t,y) 4./(1 + t.^2);
>> [t, y] =ode45(eqn,[0 1],y0);
>> y(end)
ans =

7.5.5. 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.

See also

MATLAB documentation on solving Systems of ODEs, and Higher-Order ODEs.