You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: docs/src/examples/perturbation.md
+19-3Lines changed: 19 additions & 3 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -5,46 +5,57 @@ In the [Mixed Symbolic-Numeric Perturbation Theory tutorial](https://symbolics.j
5
5
## Free fall in a varying gravitational field
6
6
7
7
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 $ẋ(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
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 $ẍ = -g$ with constant gravitational acceleration $g = GM/R^2$ and solution $x(t) = x(0) + ẋ(0) t - \frac{1}{2} g t^2$. This is a good setup for perturbation theory.
12
14
13
15
To make the problem dimensionless, we redefine $x \leftarrow x / R$ and $t \leftarrow t / \sqrt{R^3/GM}$. Then the ODE becomes
16
+
14
17
```@example perturbation
15
18
using ModelingToolkit
16
19
using ModelingToolkit: t_nounits as t, D_nounits as D
17
20
@variables ϵ x(t)
18
-
eq = D(D(x)) ~ -(1 + ϵ*x)^(-2)
21
+
eq = D(D(x)) ~ -(1 + ϵ * x)^(-2)
19
22
```
20
23
21
24
Next, expand $x(t)$ in a series up to second order in $ϵ$:
25
+
22
26
```@example perturbation
23
27
using Symbolics
24
28
@variables y(t)[0:2] # coefficients
25
29
x_series = series(y, ϵ)
26
30
```
27
31
28
32
Insert this into the equation and collect perturbed equations to each order:
33
+
29
34
```@example perturbation
30
35
eq_pert = substitute(eq, x => x_series)
31
36
eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
32
37
```
38
+
33
39
!!! note
40
+
34
41
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.
35
-
These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities:
42
+
These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities:
43
+
36
44
```@example perturbation
37
45
@mtkbuild sys = ODESystem(eqs_pert, t)
38
46
```
39
47
40
48
To solve the `ODESystem`, we generate an `ODEProblem` with initial conditions $x(0) = 0$, and $ẋ(0) = 1$, and solve it:
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 $ϵ$:
58
+
48
59
```@example perturbation
49
60
using Plots
50
61
p = plot()
@@ -53,6 +64,7 @@ for ϵᵢ in 0.0:0.1:1.0
53
64
end
54
65
p
55
66
```
67
+
56
68
This makes sense: for larger $ϵ$, gravity weakens with altitude, and the trajectory goes higher for a fixed initial velocity.
57
69
58
70
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 $ϵ$.
@@ -62,19 +74,23 @@ An advantage of the perturbative method is that we run the ODE solver only once
62
74
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.
63
75
64
76
The goal is to solve the ODE
77
+
65
78
```@example perturbation
66
-
eq = D(D(x)) + 2*ϵ*D(x) + x ~ 0
79
+
eq = D(D(x)) + 2 * ϵ * D(x) + x ~ 0
67
80
```
81
+
68
82
with initial conditions $x(0) = 0$ and $ẋ(0) = 1$. With $ϵ = 0$, the problem reduces to the simple linear harmonic oscillator with the exact solution $x(t) = \sin(t)$.
69
83
70
84
We follow the same steps as in the previous example to construct the `ODESystem`:
85
+
71
86
```@example perturbation
72
87
eq_pert = substitute(eq, x => x_series)
73
88
eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
74
89
@mtkbuild sys = ODESystem(eqs_pert, t)
75
90
```
76
91
77
92
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:
0 commit comments