Symbolic-Numeric Perturbation Theory for ODEs

In the Mixed Symbolic-Numeric Perturbation Theory tutorial, we discussed how to solve algebraic equations using Symbolics.jl. Here we extend the method to differential equations. The procedure is similar, but the Taylor series coefficients now become functions of an independent variable (usually time).

Free fall in a varying gravitational field

Our first ODE example is a well-known physics problem: what is the altitude $x (t)$ of an object (say, a ball or a rocket) thrown vertically with initial velocity $x ̇ (0)$ from the surface of a planet with mass $M$ and radius $R$ ? According to Newton's second law and law of gravity, it is the solution of the ODE

$x ̈ = - \frac{G M}{(R + x)^{2}} = - \frac{G M}{R^{2}} \frac{1}{{(1 + ϵ \frac{x}{R})}^{2}} .$

In the last equality, we introduced a perturbative expansion parameter $ϵ$ . When $ϵ = 1$ , we recover the original problem. When $ϵ = 0$ , the problem reduces to the trivial problem $x ̈ = - g$ with constant gravitational acceleration $g = G M / R^{2}$ and solution $x (t) = x (0) + x ̇ (0) t - \frac{1}{2} g t^{2}$ . This is a good setup for perturbation theory.

To make the problem dimensionless, we redefine $x \leftarrow x / R$ and $t \leftarrow t / \sqrt{R^{3} / G M}$ . Then the ODE becomes

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables ϵ x(t)
eq = D(D(x)) ~ -(1 + ϵ * x)^(-2)

$\frac{d}{d t} \frac{d x (t)}{d t} = - {(\frac{1}{1 + x (t) ϵ})}^{2}$

Next, expand $x (t)$ in a series up to second order in $ϵ$ :

using Symbolics
@variables y(t)[0:2] # coefficients
x_series = series(y, ϵ)

$y_0 (t) + y_1 (t) ϵ + ϵ^{2} y_2 (t)$

Insert this into the equation and collect perturbed equations to each order:

eq_pert = substitute(eq, x => x_series)
eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)

$\begin{aligned} \frac{d}{d t} \frac{d y_0 (t)}{d t} & = - 1 \\ \frac{d}{d t} \frac{d y_1 (t)}{d t} & = 2 y_0 (t) \\ \frac{d}{d t} \frac{d y_2 (t)}{d t} & = \frac{1}{2} (4 y_1 (t) - 6 y_0 {(t)}^{2}) \end{aligned}$

Note

The 0-th order equation can be solved analytically, but ModelingToolkit does currently not feature automatic analytical solution of ODEs, so we proceed with solving it numerically.

These are the ODEs we want to solve. Now construct an ODESystem, which automatically inserts dummy derivatives for the velocities:

@mtkbuild sys = ODESystem(eqs_pert, t)

$\begin{aligned} \frac{d y ˍ t_0 (t)}{d t} & = - 1 \\ \frac{d y ˍ t_1 (t)}{d t} & = 2 y_0 (t) \\ \frac{d y ˍ t_2 (t)}{d t} & = 0.5 (4 y_1 (t) - 6 y_0 {(t)}^{2}) \\ \frac{d y_0 (t)}{d t} & = y ˍ t_0 (t) \\ \frac{d y_1 (t)}{d t} & = y ˍ t_1 (t) \\ \frac{d y_2 (t)}{d t} & = y ˍ t_2 (t) \end{aligned}$

To solve the ODESystem, we generate an ODEProblem with initial conditions $x (0) = 0$ , and $x ̇ (0) = 1$ , and solve it:

using OrdinaryDiffEq
u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
prob = ODEProblem(sys, u0, (0.0, 3.0))
sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 17-element Vector{Float64}:
 0.0
 0.0009990005004983772
 0.010989005505482149
 0.050323201520099475
 0.10988625821422165
 0.17976131889392957
 0.2685822002061843
 0.37704291021657915
 0.5150699918763888
 0.6895893855495662
 0.9121374549979837
 1.1969237822790968
 1.5496867745920098
 1.9254996658947618
 2.3244234888380926
 2.747063109192439
 3.0
u: 17-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 [0.9990009994995016, 9.976696651635089e-7, -9.96091672587352e-10, 0.0009985014994983792, 3.3225183183212947e-10, -2.488196097577267e-13]
 [0.9890109944945178, 0.00012031590433746151, -1.313675024660254e-6, 0.010928626384482389, 4.4112244930411036e-7, -3.616313645853468e-9]
 [0.9496767984799005, 0.0024899447065507025, -0.00012162013827075771, 0.049056989214483204, 4.1945473480873235e-5, -1.5446236873125787e-6]
 [0.8901137417857783, 0.01163269793066302, -0.0011961578797763368, 0.10384876334206036, 0.00043014136554910947, -3.356778897522297e-5]
 [0.8202386811060703, 0.030377854788493475, -0.004886057770148157, 0.16360425300868706, 0.0018492600559848445, -0.00022766884565584688]
 [0.7314177997938156, 0.06567821408037469, -0.014860759516862308, 0.23251400107238687, 0.0060245458576349475, -0.0010561563724513826]
 [0.6229570897834208, 0.12429437899755692, -0.036472231061079205, 0.3059622321442855, 0.016182822881957857, -0.0037432549829311493]
 [0.4849300081236112, 0.2197482387464324, -0.0787754252259981, 0.3824214436106172, 0.03968364533277146, -0.011519993739428212]
 [0.3104106144504339, 0.3662258978886398, -0.1492237744354559, 0.4518226252182521, 0.0904632786678382, -0.03123009689695499]
 [0.0878625450020164, 0.5790302162082741, -0.2401187555898396, 0.49614008659288455, 0.19527991709526749, -0.07489557108797736]
 [-0.19692378227909657, 0.8610449480683734, -0.2837355009128571, 0.48061051198644794, 0.4005466921178069, -0.15257432818111172]
 [-0.5496867745920095, 1.160989804660918, -0.07344287192068476, 0.34892222491931746, 0.7599274601009123, -0.2264891132383536]
 [-0.9254996658947615, 1.3279208666142521, 0.6091185233659987, 0.07172518421434361, 1.2341348204367515, -0.14128168530832333]
 [-1.324423488838092, 1.2167008109268078, 1.7605042472862638, -0.37704878889303023, 1.753592922085515, 0.3226927609070924]
 [-1.747063109192438, 0.6362505847444941, 2.791184265640095, -1.0261147537505726, 2.164481412673811, 1.3123578427272624]
 [-1.9999999999999991, 1.426739902911473e-14, 2.700000000000024, -1.4999999999999962, 2.2500000000000058, 2.02499470614646]

This is the solution for the coefficients in the series for $x (t)$ and their derivatives. Finally, we calculate the solution to the original problem by summing the series for different $ϵ$ :

using Plots
p = plot()
for ϵᵢ in 0.0:0.1:1.0
    plot!(p, sol, idxs = substitute(x_series, ϵ => ϵᵢ), label = "ϵ = $ϵᵢ")
end
p

This makes sense: for larger $ϵ$ , gravity weakens with altitude, and the trajectory goes higher for a fixed initial velocity.

An advantage of the perturbative method is that we run the ODE solver only once and calculate trajectories for several $ϵ$ for free. Had we solved the full unperturbed ODE directly, we would need to do repeat it for every $ϵ$ .

Weakly nonlinear oscillator

Our second example applies perturbation theory to nonlinear oscillators – a very important class of problems. As we will see, perturbation theory has difficulty providing a good solution to this problem, but the process is nevertheless instructive. This example closely follows chapter 7.6 of Nonlinear Dynamics and Chaos by Steven Strogatz.

The goal is to solve the ODE

eq = D(D(x)) + 2 * ϵ * D(x) + x ~ 0

$x (t) + \frac{d}{d t} \frac{d x (t)}{d t} + 2 \frac{d x (t)}{d t} ϵ = 0$

with initial conditions $x (0) = 0$ and $x ̇ (0) = 1$ . With $ϵ = 0$ , the problem reduces to the simple linear harmonic oscillator with the exact solution $x (t) = \sin (t)$ .

We follow the same steps as in the previous example to construct the ODESystem:

eq_pert = substitute(eq, x => x_series)
eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
@mtkbuild sys = ODESystem(eqs_pert, t)

$\begin{aligned} \frac{d y_1 (t)}{d t} & = y ˍ t_1 (t) \\ \frac{d y ˍ t_2 (t)}{d t} & = - \frac{1}{2} (2 y_2 (t) + 4 y ˍ t_1 (t)) \\ \frac{d y ˍ t_0 (t)}{d t} & = - y_0 (t) \\ \frac{d y_0 (t)}{d t} & = y ˍ t_0 (t) \\ \frac{d y ˍ t_1 (t)}{d t} & = - y_1 (t) - 2 y ˍ t_0 (t) \\ \frac{d y_2 (t)}{d t} & = y ˍ t_2 (t) \end{aligned}$

We solve and plot it as in the previous example, and compare the solution with $ϵ = 0.1$ to the exact solution $x (t, ϵ) = e^{- ϵ t} \sin (\sqrt{(1 - ϵ^{2})} t) / \sqrt{1 - ϵ^{2}}$ of the unperturbed equation:

u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
prob = ODEProblem(sys, u0, (0.0, 50.0))
sol = solve(prob)
plot(sol, idxs = substitute(x_series, ϵ => 0.1); label = "Perturbative (ϵ=0.1)")

x_exact(t, ϵ) = exp(-ϵ * t) * sin(√(1 - ϵ^2) * t) / √(1 - ϵ^2)
plot!(sol.t, x_exact.(sol.t, 0.1); label = "Exact (ϵ=0.1)")

This is similar to Figure 7.6.2 in Nonlinear Dynamics and Chaos. The two curves fit well for the first couple of cycles, but then the perturbative solution diverges from the exact solution. The main reason is that the problem has two or more time-scales that introduce secular terms in the solution. One solution is to explicitly account for the two time scales and use an analytic method called two-timing, but this is outside the scope of this example.