Ordinary differential equations

Solving the ODE initial value problem (odefun)

mpmath.odefun(F, x0, y0, tol=None, degree=None, method='taylor', verbose=False)

Returns a function \(y(x) = [y_0(x), y_1(x), \ldots, y_n(x)]\) that is a numerical solution of the \(n+1\)-dimensional first-order ordinary differential equation (ODE) system

\[ \begin{align}\begin{aligned}y_0'(x) = F_0(x, [y_0(x), y_1(x), \ldots, y_n(x)])\\y_1'(x) = F_1(x, [y_0(x), y_1(x), \ldots, y_n(x)])\\\vdots\\y_n'(x) = F_n(x, [y_0(x), y_1(x), \ldots, y_n(x)])\end{aligned}\end{align} \]

The derivatives are specified by the vector-valued function F that evaluates \([y_0', \ldots, y_n'] = F(x, [y_0, \ldots, y_n])\). The initial point \(x_0\) is specified by the scalar argument x0, and the initial value \(y(x_0) = [y_0(x_0), \ldots, y_n(x_0)]\) is specified by the vector argument y0.

For convenience, if the system is one-dimensional, you may optionally provide just a scalar value for y0. In this case, F should accept a scalar y argument and return a scalar. The solution function y will return scalar values instead of length-1 vectors.

Evaluation of the solution function \(y(x)\) is permitted for any \(x \ge x_0\).

A high-order ODE can be solved by transforming it into first-order vector form. This transformation is described in standard texts on ODEs. Examples will also be given below.

Options, speed and accuracy

By default, odefun() uses a high-order Taylor series method. For reasonably well-behaved problems, the solution will be fully accurate to within the working precision. Note that F must be possible to evaluate to very high precision for the generation of Taylor series to work.

To get a faster but less accurate solution, you can set a large value for tol (which defaults roughly to eps). If you just want to plot the solution or perform a basic simulation, tol = 0.01 is likely sufficient.

The degree argument controls the degree of the solver (with method=’taylor’, this is the degree of the Taylor series expansion). A higher degree means that a longer step can be taken before a new local solution must be generated from F, meaning that fewer steps are required to get from \(x_0\) to a given \(x_1\). On the other hand, a higher degree also means that each local solution becomes more expensive (i.e., more evaluations of F are required per step, and at higher precision).

The optimal setting therefore involves a tradeoff. Generally, decreasing the degree for Taylor series is likely to give faster solution at low precision, while increasing is likely to be better at higher precision.

The function object returned by odefun() caches the solutions at all step points and uses polynomial interpolation between step points. Therefore, once \(y(x_1)\) has been evaluated for some \(x_1\), \(y(x)\) can be evaluated very quickly for any \(x_0 \le x \le x_1\). and continuing the evaluation up to \(x_2 > x_1\) is also fast.

Examples of first-order ODEs

We will solve the standard test problem \(y'(x) = y(x), y(0) = 1\) which has explicit solution \(y(x) = \exp(x)\):

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> f = odefun(lambda x, y: y, 0, 1)
>>> for x in [0, 1, 2.5]:
...     print((f(x), exp(x)))
...
(1.0, 1.0)
(2.71828182845905, 2.71828182845905)
(12.1824939607035, 12.1824939607035)

The solution with high precision:

>>> mp.dps = 50
>>> f = odefun(lambda x, y: y, 0, 1)
>>> f(1)
2.7182818284590452353602874713526624977572470937
>>> exp(1)
2.7182818284590452353602874713526624977572470937

Using the more general vectorized form, the test problem can be input as (note that f returns a 1-element vector):

>>> mp.dps = 15
>>> f = odefun(lambda x, y: [y[0]], 0, [1])
>>> f(1)
[2.71828182845905]

odefun() can solve nonlinear ODEs, which are generally impossible (and at best difficult) to solve analytically. As an example of a nonlinear ODE, we will solve \(y'(x) = x \sin(y(x))\) for \(y(0) = \pi/2\). An exact solution happens to be known for this problem, and is given by \(y(x) = 2 \tan^{-1}\left(\exp\left(x^2/2\right)\right)\):

>>> f = odefun(lambda x, y: x*sin(y), 0, pi/2)
>>> for x in [2, 5, 10]:
...     print((f(x), 2*atan(exp(mpf(x)**2/2))))
...
(2.87255666284091, 2.87255666284091)
(3.14158520028345, 3.14158520028345)
(3.14159265358979, 3.14159265358979)

If \(F\) is independent of \(y\), an ODE can be solved using direct integration. We can therefore obtain a reference solution with quad():

>>> f = lambda x: (1+x**2)/(1+x**3)
>>> g = odefun(lambda x, y: f(x), pi, 0)
>>> g(2*pi)
0.72128263801696
>>> quad(f, [pi, 2*pi])
0.72128263801696

Examples of second-order ODEs

We will solve the harmonic oscillator equation \(y''(x) + y(x) = 0\). To do this, we introduce the helper functions \(y_0 = y, y_1 = y_0'\) whereby the original equation can be written as \(y_1' + y_0' = 0\). Put together, we get the first-order, two-dimensional vector ODE

\[\begin{split}\begin{cases} y_0' = y_1 \\ y_1' = -y_0 \end{cases}\end{split}\]

To get a well-defined IVP, we need two initial values. With \(y(0) = y_0(0) = 1\) and \(-y'(0) = y_1(0) = 0\), the problem will of course be solved by \(y(x) = y_0(x) = \cos(x)\) and \(-y'(x) = y_1(x) = \sin(x)\). We check this:

>>> f = odefun(lambda x, y: [-y[1], y[0]], 0, [1, 0])
>>> for x in [0, 1, 2.5, 10]:
...     nprint(f(x), 15)
...     nprint([cos(x), sin(x)], 15)
...     print("---")
...
[1.0, 0.0]
[1.0, 0.0]
---
[0.54030230586814, 0.841470984807897]
[0.54030230586814, 0.841470984807897]
---
[-0.801143615546934, 0.598472144103957]
[-0.801143615546934, 0.598472144103957]
---
[-0.839071529076452, -0.54402111088937]
[-0.839071529076452, -0.54402111088937]
---

Note that we get both the sine and the cosine solutions simultaneously.

TODO

  • Better automatic choice of degree and step size

  • Make determination of Taylor series convergence radius more robust

  • Allow solution for \(x < x_0\)

  • Allow solution for complex \(x\)

  • Test for difficult (ill-conditioned) problems

  • Implement Runge-Kutta and other algorithms