Skip to content

Commit 7ca38cb

Browse files
refactor: format
1 parent 5a6da3c commit 7ca38cb

File tree

2 files changed

+21
-5
lines changed

2 files changed

+21
-5
lines changed

docs/src/examples/perturbation.md

Lines changed: 19 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -5,46 +5,57 @@ In the [Mixed Symbolic-Numeric Perturbation Theory tutorial](https://symbolics.j
55
## Free fall in a varying gravitational field
66

77
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
8+
89
```math
910
ẍ = -\frac{GM}{(R+x)^2} = -\frac{GM}{R^2} \frac{1}{\left(1+ϵ\frac{x}{R}\right)^2}.
1011
```
12+
1113
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.
1214

1315
To make the problem dimensionless, we redefine $x \leftarrow x / R$ and $t \leftarrow t / \sqrt{R^3/GM}$. Then the ODE becomes
16+
1417
```@example perturbation
1518
using ModelingToolkit
1619
using ModelingToolkit: t_nounits as t, D_nounits as D
1720
@variables ϵ x(t)
18-
eq = D(D(x)) ~ -(1 + ϵ*x)^(-2)
21+
eq = D(D(x)) ~ -(1 + ϵ * x)^(-2)
1922
```
2023

2124
Next, expand $x(t)$ in a series up to second order in $ϵ$:
25+
2226
```@example perturbation
2327
using Symbolics
2428
@variables y(t)[0:2] # coefficients
2529
x_series = series(y, ϵ)
2630
```
2731

2832
Insert this into the equation and collect perturbed equations to each order:
33+
2934
```@example perturbation
3035
eq_pert = substitute(eq, x => x_series)
3136
eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
3237
```
38+
3339
!!! note
40+
3441
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+
3644
```@example perturbation
3745
@mtkbuild sys = ODESystem(eqs_pert, t)
3846
```
3947

4048
To solve the `ODESystem`, we generate an `ODEProblem` with initial conditions $x(0) = 0$, and $ẋ(0) = 1$, and solve it:
49+
4150
```@example perturbation
4251
using DifferentialEquations
4352
u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
4453
prob = ODEProblem(sys, u0, (0.0, 3.0))
4554
sol = solve(prob)
4655
```
56+
4757
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+
4859
```@example perturbation
4960
using Plots
5061
p = plot()
@@ -53,6 +64,7 @@ for ϵᵢ in 0.0:0.1:1.0
5364
end
5465
p
5566
```
67+
5668
This makes sense: for larger $ϵ$, gravity weakens with altitude, and the trajectory goes higher for a fixed initial velocity.
5769

5870
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
6274
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.
6375

6476
The goal is to solve the ODE
77+
6578
```@example perturbation
66-
eq = D(D(x)) + 2*ϵ*D(x) + x ~ 0
79+
eq = D(D(x)) + 2 * ϵ * D(x) + x ~ 0
6780
```
81+
6882
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)$.
6983

7084
We follow the same steps as in the previous example to construct the `ODESystem`:
85+
7186
```@example perturbation
7287
eq_pert = substitute(eq, x => x_series)
7388
eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
7489
@mtkbuild sys = ODESystem(eqs_pert, t)
7590
```
7691

7792
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:
93+
7894
```@example perturbation
7995
u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
8096
prob = ODEProblem(sys, u0, (0.0, 50.0))

src/systems/abstractsystem.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3001,8 +3001,8 @@ By default, the resulting system inherits `sys`'s name and description.
30013001
See also [`compose`](@ref).
30023002
"""
30033003
function extend(sys::AbstractSystem, basesys::AbstractSystem;
3004-
name::Symbol = nameof(sys), description = description(sys),
3005-
gui_metadata = get_gui_metadata(sys))
3004+
name::Symbol = nameof(sys), description = description(sys),
3005+
gui_metadata = get_gui_metadata(sys))
30063006
T = SciMLBase.parameterless_type(basesys)
30073007
ivs = independent_variables(basesys)
30083008
if !(sys isa T)

0 commit comments

Comments
 (0)