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
+17-17Lines changed: 17 additions & 17 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -4,7 +4,7 @@
4
4
5
5
In the previous tutorial, [Mixed Symbolic-Numeric Perturbation Theory](https://symbolics.juliasymbolics.org/stable/examples/perturbation), we discussed how to solve algebraic equations using **Symbolics.jl**. Here, our goal is to extend the method to differential equations. First, we import the following helper functions that were introduced in [Mixed Symbolic/Numerical Methods for Perturbation Theory - Algebraic Equations](@ref perturb_alg):
6
6
7
-
```@example perturb
7
+
```julia
8
8
using Symbolics, SymbolicUtils
9
9
10
10
def_taylor(x, ps) =sum([a*x^i for (i,a) inenumerate(ps)])
@@ -44,60 +44,60 @@ As the first ODE example, we have chosen a simple and well-behaved problem, whic
44
44
45
45
with the initial conditions $x(0) = 0$, and $\dot{x}(0) = 1$. Note that for $\epsilon = 0$, this equation transforms back to the standard one. Let's start with defining the variables
46
46
47
-
```@example perturb
47
+
```julia
48
48
n =3
49
49
@variables ϵ t y[1:n](t) ∂∂y[1:n](t)
50
50
```
51
51
52
52
Next, we define $x$.
53
53
54
-
```@example perturb
54
+
```julia
55
55
x =def_taylor(ϵ, y[3:end], y[2])
56
56
```
57
57
58
58
We need the second derivative of `x`. It may seem that we can do this using `Differential(t)`; however, this operation needs to wait for a few steps because we need to manipulate the differentials as separate variables. Instead, we define dummy variables `∂∂y` as the placeholder for the second derivatives and define
59
59
60
-
```@example perturb
60
+
```julia
61
61
∂∂x =def_taylor(ϵ, ∂∂y[3:end], ∂∂y[2])
62
62
```
63
63
64
64
as the second derivative of `x`. After rearrangement, our governing equation is $\ddot{x}(t)(1 + \epsilon x(t))^{-2} + 1 = 0$, or
65
65
66
-
```@example perturb
66
+
```julia
67
67
eq = ∂∂x * (1+ ϵ*x)^2+1
68
68
```
69
69
70
70
The next two steps are the same as the ones for algebraic equations (note that we pass `1:n` to `collect_powers` because the zeroth order term is needed here)
71
71
72
-
```@example perturb
72
+
```julia
73
73
eqs =collect_powers(eq, ϵ, 1:n)
74
74
```
75
75
76
76
and,
77
77
78
-
```@example perturb
78
+
```julia
79
79
vals =solve_coef(eqs, ∂∂y)
80
80
```
81
81
82
82
Our system of ODEs is forming. Now is the time to convert `∂∂`s to the correct **Symbolics.jl** form by substitution:
83
83
84
-
```@example perturb
84
+
```julia
85
85
D =Differential(t)
86
86
subs =Dict(∂∂y[i] =>D(D(y[i])) for i ineachindex(y))
87
87
eqs = [substitute(first(v), subs) ~substitute(last(v), subs) for v in vals]
88
88
```
89
89
90
90
We are nearly there! From this point on, the rest is standard ODE solving procedures. Potentially we can use a symbolic ODE solver to find a closed form solution to this problem. However, **Symbolics.jl** currently does not support this functionality. Instead, we solve the problem numerically. We form an `ODESystem`, lower the order (convert second derivatives to first), generate an `ODEProblem` (after passing the correct initial conditions), and, finally, solve it.
91
91
92
-
```@example perturb
92
+
```julia
93
93
using ModelingToolkit, DifferentialEquations
94
94
95
95
@named sys =ODESystem(eqs, t)
96
96
sys =structural_simplify(sys)
97
97
states(sys)
98
98
```
99
99
100
-
```@example perturb
100
+
```julia
101
101
# the initial conditions
102
102
# everything is zero except the initial velocity
103
103
u0 =zeros(2n+2)
@@ -109,13 +109,13 @@ sol = solve(prob; dtmax=0.01)
109
109
110
110
Finally, we calculate the solution to the problem as a function of `ϵ` by substituting the solution to the ODE system back into the defining equation for `x`. Note that `𝜀` is a number, compared to `ϵ`, which is a symbolic variable.
111
111
112
-
```@example perturb
112
+
```julia
113
113
X = 𝜀 ->sum([𝜀^(i-1) * sol[y[i]] for i ineachindex(y)])
114
114
```
115
115
116
116
Using `X`, we can plot the trajectory for a range of $𝜀$s.
117
117
118
-
```@example perturb
118
+
```julia
119
119
using Plots
120
120
121
121
plot(sol.t, hcat([X(𝜀) for 𝜀 =0.0:0.1:0.5]...))
@@ -129,7 +129,7 @@ For the next example, we have chosen a simple example from a very important clas
129
129
130
130
The goal is to solve $\ddot{x} + 2\epsilon\dot{x} + x = 0$, where the dot signifies time-derivatives and the initial conditions are $x(0) = 0$ and $\dot{x}(0) = 1$. If $\epsilon = 0$, the problem reduces to the simple linear harmonic oscillator with the exact solution $x(t) = \sin(t)$. We follow the same steps as the previous example.
131
131
132
-
```@example perturb
132
+
```julia
133
133
n =3
134
134
@variables ϵ t y[1:n](t) ∂y[1:n] ∂∂y[1:n]
135
135
x =def_taylor(ϵ, y[3:end], y[2])
@@ -139,15 +139,15 @@ x = def_taylor(ϵ, y[3:end], y[2])
139
139
140
140
This time we also need the first derivative terms. Continuing,
141
141
142
-
```@example perturb
142
+
```julia
143
143
eq = ∂∂x +2*ϵ*∂x + x
144
144
eqs =collect_powers(eq, ϵ, 0:n)
145
145
vals =solve_coef(eqs, ∂∂y)
146
146
```
147
147
148
148
Next, we need to replace `∂`s and `∂∂`s with their **Symbolics.jl** counterparts:
149
149
150
-
```@example perturb
150
+
```julia
151
151
D =Differential(t)
152
152
subs1 =Dict(∂y[i] =>D(y[i]) for i ineachindex(y))
153
153
subs2 =Dict(∂∂y[i] =>D(D(y[i])) for i ineachindex(y))
@@ -157,12 +157,12 @@ eqs = [substitute(first(v), subs) ~ substitute(last(v), subs) for v in vals]
157
157
158
158
We continue with converting 'eqs' to an `ODEProblem`, solving it, and finally plot the results against the exact solution to the original problem, which is $x(t, \epsilon) = (1 - \epsilon)^{-1/2} e^{-\epsilon t} \sin((1- \epsilon^2)^{1/2}t)$,
0 commit comments