python22.sci2u.dk Active sci2u
Loading...
Login
Prefer latest booklet Prefer latest booklet Login
Link to booklet is copied to clipboard!
Saved to booklet!
Removed from booklet!
Started capturing to booklet.
Stopped capturing to booklet.
Preferring latest booklet.
No longer prefers latest booklet.
Entered edit mode.
Exited edit mode.
Edit mode

Click on the booklet title or a chapter title to edit (max. 30 characters)


How To's

How to get python How to plot points and lines How to use loops How to draw general shapes How to draw curves How to use functions How to manipulate Polygons How to create animations How to do mathematical analysis How to fit How to go from SymPy to NumPy How to solve a single differential equation
12.1 A single ODE
12.2 ODE: More datapoints
12.3 ODE: Backwards integration
12.4 ODE: Accuracy
How to do coupled and 2nd order ODEs How to combine vectors and ODEs How to do ODE events Add Chapter..

How to solve a single differential equation

12.1 A single ODE

Here is an example of a first order ordinary differential equation (ODE):
The ODE is a so-called initial value problem (IVP) whose solution is well-defined once the initial value condition is specified. That may e.g. be that the point:
is part of the solution.
Assume we want to solve the ODE for -values fulfilling , we may start by declaring a list (or tuple or NumPy array) with the initial and final -values:
tinit = 0
tfinal = 3
trange = [tinit,tfinal]
and a list (or NumPy array, but not a tuple) with the -value dictated by the initial value condition:
yinit = [1]
At this point, it may seem unnecessary to declare this value as a list, but that will turn out logical once we are to solve several differential equations at the same time.
As a final preparation, let us define a function that given and evaluates :
import numpy as np
dydt = lambda t, y: 0.9 * y**2 * np.cos(t)
We are now ready to import solve_ivp from scipy.integrate and call it like this:
from scipy.integrate import solve_ivp

mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
ys = mysol.y[0]
There are three arguments for solve_ivp:
  1. dydt is the function that implements the differential equation
  2. trange is a list with the starting and ending -values for which the differential equation should be solved. Importantly, the first element of trange has to conform with that of the initial values.
  3. yinit is a list with the -value dictated by the initial values.
The ts = mysol.t and ys = mysol.y[0] following the call of solve_ivp select from the output of the function, the and values for which the differential equation was actually solved.
These values may then be plotted with:
import matplotlib.pyplot as plt

plt.rc('font', size=16)
fig,ax = plt.subplots(1,1)
ax.plot(ts,ys,'o',linestyle='--')
ax.grid()
ax.set_ylim([0,11])
ax.set_xlabel('$t$')
ax.set_ylabel('$y$')

Figure 1



12.2 ODE: More datapoints

Using solve_ivp to solve a differential equation, the user may specify for which -values, the solution is sought. That is done by adding a keyword-argument t_eval with a Python list (or NumPy array) containing the values:
ts = np.linspace(tinit, tfinal, 100)
mysol = solve_ivp(dydt,
                  trange,
                  yinit,
                  t_eval=ts)
ts = mysol.t
ys = mysol.y[0]

ax.plot(ts,ys)
fig

Figure 2



12.3 ODE: Backwards integration

Assume the differential equation is to be solved with this initial value condition:
and for -values fulfilling , i.e. for values before that of the initial condition. Well, it makes no difference to solve_ivp, that the integration now takes place from large values to smaller ones. Only, it is important that the first -value given conforms with that of the initial value condition. Thus, the trange and yinit must be changed into:
tinit = 9
tfinal = 4
trange = (tinit,tfinal)
yinit = [5/4]
After which, the same call of solve_ivp can be made and the result be plotted:
mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
ys = mysol.y[0]

ax.plot(ts,ys,'o',linestyle='--')
ax.set_xlim([-1,10])
fig
Note how trange is defined so that the first value is the one conforming with the new initial values.

Figure 3



12.4 ODE: Accuracy

The differential equation above can be solve by separation of the variables and integration. It turns out to have the exact solution:
Using initial values, , and solving until one should thus get:
as the final value for . Lets inspect that:
tinit = 0
tfinal = np.pi
trange = [tinit,tfinal]
yinit = [1]
mysol = solve_ivp(dydt, trange, yinit)

ys = mysol.y[0]

# now print the final value for y, corresponding to y(tfinal) = y(t = pi)
ys[-1]
1.0107395513007404
Well, that was not exactly the case. This is due to rounding errors made by solve_ivp while solving the differential equation.
One may control the accuracy by setting e.g. a maximum step size for the independent variable, here , during the integration:
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-1)

ys = mysol.y[0]

ys[-1]
1.0000016575704107
Here the maximum step size was set to and the numerical result for came much closer to exact .
Using np.logspace:
np.logspace(0,-4,5)
array([1.e+00, 1.e-01, 1.e-02, 1.e-03, 1.e-04])
and a list comprehension we may construct a neat list of max_step-values and -values in this way:
res = [(max_step, solve_ivp(dydt,trange,yinit,max_step=max_step).y[0][-1]) \
             for max_step in np.logspace(0,-4,5)]
res
Perhaps not so easy to read. You might prefer to write it like this:
res = []
for max_step in np.logspace(0,-4,5):
  mysol = solve_ivp(dydt, trange, yinit, max_step=max_step)
  ys = mysol.y[0]
  res.append((max_step, ys[-1]))
res
They both give the exact same result, namely:
[(1.0, 1.0107395513007404),
 (0.1, 1.0000016575704107),
 (0.01, 1.0000000000079718),
 (0.001, 0.9999999999999999),
 (0.0001, 0.9999999999999947)]
i.e. practically with a sufficiently small value for max_step.
Other arguments can be set when calling solve_ivp. These include the relative tolerance, rtol, and the absolute tolerance, atol. You may find more information here.


Sci2u Assignment: 820
Delete "How To's"?
Once deleted this booklet is gone forever!
Block is found in multiple booklets!

Choose which booklet to go to:

© 2019-2022 Uniblender ApS