# 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:

Considered in terms of discrete samples of , the derivative becomes

Rearranged to find given the previous value, , and the slope of the function, , we have

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, and initial conditions, .

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 , 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 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

Then `[t,y] = ode45(f,tspan,y0)`

, where `tspan = [t0 tf]`

,
integrates the system of differential equations from
to with initial conditions . 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 .

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

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

Where, .

In the distant past, (), the population size is modeled as zero. The steady state () population size is , which is the sustainable population size.

The half way point, can be modeled as .

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 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 , then

For example, consider

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
```

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