@@ -19,7 +19,7 @@ x' = f(x,y,t)\\
19
19
```
20
20
21
21
where `` x `` are the differential variables and `` y `` are the algebraic variables.
22
- An initial condition `` u0 = [x(t_0) y(t_0)] `` is said to be consistent if
22
+ An initial condition `` u0 = [x(t_0) y(t_0)] `` is said to be consistent if
23
23
`` g(x(t_0),y(t_0),t_0) = 0 `` .
24
24
25
25
For ODEs, this is trivially satisfied. However, for more complicated systems it may
@@ -64,13 +64,13 @@ This solves via:
64
64
65
65
``` @example init
66
66
sol = solve(prob, Rodas5P())
67
- plot(sol, idxs = (x,y))
67
+ plot(sol, idxs = (x, y))
68
68
```
69
69
70
70
and we can check it satisfies our conditions via:
71
71
72
72
``` @example init
73
- conditions = getfield.(equations(pend)[3:end],:rhs)
73
+ conditions = getfield.(equations(pend)[3:end], :rhs)
74
74
```
75
75
76
76
``` @example init
@@ -93,48 +93,50 @@ We can similarly choose `λ = 0` and solve for `y` to start the system:
93
93
``` @example init
94
94
prob = ODEProblem(pend, [x => 1, λ => 0], (0.0, 1.5), [g => 1], guesses = [y => 1])
95
95
sol = solve(prob, Rodas5P())
96
- plot(sol, idxs = (x,y))
96
+ plot(sol, idxs = (x, y))
97
97
```
98
98
99
99
or choose to satisfy derivative conditions:
100
100
101
101
``` @example init
102
- prob = ODEProblem(pend, [x => 1, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1])
102
+ prob = ODEProblem(
103
+ pend, [x => 1, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1])
103
104
sol = solve(prob, Rodas5P())
104
- plot(sol, idxs = (x,y))
105
+ plot(sol, idxs = (x, y))
105
106
```
106
107
107
- Notice that since a derivative condition is given, we are required to give a
108
+ Notice that since a derivative condition is given, we are required to give a
108
109
guess for ` y ` .
109
110
110
111
We can also directly give equations to be satisfied at the initial point by using
111
112
the ` initialization_eqs ` keyword argument, for example:
112
113
113
114
``` @example init
114
115
prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1],
115
- initialization_eqs = [y ~ 1])
116
+ initialization_eqs = [y ~ 1])
116
117
sol = solve(prob, Rodas5P())
117
- plot(sol, idxs = (x,y))
118
+ plot(sol, idxs = (x, y))
118
119
```
119
120
120
121
Additionally, note that the initial conditions are allowed to be functions of other
121
122
variables and parameters:
122
123
123
124
``` @example init
124
- prob = ODEProblem(pend, [x => 1, D(y) => g], (0.0, 3.0), [g => 1], guesses = [λ => 0, y => 1])
125
+ prob = ODEProblem(
126
+ pend, [x => 1, D(y) => g], (0.0, 3.0), [g => 1], guesses = [λ => 0, y => 1])
125
127
sol = solve(prob, Rodas5P())
126
- plot(sol, idxs = (x,y))
128
+ plot(sol, idxs = (x, y))
127
129
```
128
130
129
131
## Determinability: Underdetermined and Overdetermined Systems
130
132
131
133
For this system we have 3 conditions to satisfy:
132
134
133
135
``` @example init
134
- conditions = getfield.(equations(pend)[3:end],:rhs)
136
+ conditions = getfield.(equations(pend)[3:end], :rhs)
135
137
```
136
138
137
- when we initialize with
139
+ when we initialize with
138
140
139
141
``` @example init
140
142
prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = [y => 0, λ => 1])
@@ -155,21 +157,22 @@ and thus the solution is not necessarily unique. It can still be solved:
155
157
156
158
``` @example init
157
159
sol = solve(prob, Rodas5P())
158
- plot(sol, idxs = (x,y))
160
+ plot(sol, idxs = (x, y))
159
161
```
160
162
161
163
and the found initial condition satisfies all constraints which were given. In the opposite
162
164
direction, we may have an overdetermined system:
163
165
164
166
``` @example init
165
- prob = ODEProblem(pend, [x => 1, y => 0.0, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 1])
167
+ prob = ODEProblem(
168
+ pend, [x => 1, y => 0.0, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 1])
166
169
```
167
170
168
171
Can that be solved?
169
172
170
173
``` @example init
171
174
sol = solve(prob, Rodas5P())
172
- plot(sol, idxs = (x,y))
175
+ plot(sol, idxs = (x, y))
173
176
```
174
177
175
178
Indeed since we saw ` D(y) = 0 ` at the initial point above, it turns out that this solution
@@ -178,7 +181,8 @@ aren't that lucky. If the set of initial conditions cannot be satisfied, then yo
178
181
a ` SciMLBase.ReturnCode.InitialFailure ` :
179
182
180
183
``` @example init
181
- prob = ODEProblem(pend, [x => 1, y => 0.0, D(y) => 2.0, λ => 1], (0.0, 1.5), [g => 1], guesses = [λ => 1])
184
+ prob = ODEProblem(
185
+ pend, [x => 1, y => 0.0, D(y) => 2.0, λ => 1], (0.0, 1.5), [g => 1], guesses = [λ => 1])
182
186
sol = solve(prob, Rodas5P())
183
187
```
184
188
@@ -216,18 +220,19 @@ with observables, those observables are too treated as initial equations. We can
216
220
resulting simplified system via the command:
217
221
218
222
``` @example init
219
- isys = structural_simplify(isys; fully_determined= false)
223
+ isys = structural_simplify(isys; fully_determined = false)
220
224
```
221
225
222
226
Note ` fully_determined=false ` allows for the simplification to occur when the number of equations
223
227
does not match the number of unknowns, which we can use to investigate our overdetermined system:
224
228
225
229
``` @example init
226
- isys = ModelingToolkit.generate_initializesystem(pend, u0map = [x => 1, y => 0.0, D(y) => 2.0, λ => 1], guesses = [λ => 1])
230
+ isys = ModelingToolkit.generate_initializesystem(
231
+ pend, u0map = [x => 1, y => 0.0, D(y) => 2.0, λ => 1], guesses = [λ => 1])
227
232
```
228
233
229
234
``` @example init
230
- isys = structural_simplify(isys; fully_determined= false)
235
+ isys = structural_simplify(isys; fully_determined = false)
231
236
```
232
237
233
238
``` @example init
@@ -252,8 +257,8 @@ constructor which acts just like an `ODEProblem` or `NonlinearProblem` construct
252
257
creates the special initialization system for a given ` sys ` . This is done as follows:
253
258
254
259
``` @example init
255
- iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
256
- [x => 1, y => 0.0, D(y) => 2.0, λ => 1], [g => 1], guesses = [λ => 1])
260
+ iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
261
+ [x => 1, y => 0.0, D(y) => 2.0, λ => 1], [g => 1], guesses = [λ => 1])
257
262
```
258
263
259
264
We can see that because the system is overdetermined we receive a NonlinearLeastSquaresProblem,
@@ -266,6 +271,7 @@ sol = solve(iprob)
266
271
```
267
272
268
273
!!! note
274
+
269
275
For more information on solving NonlinearProblems and NonlinearLeastSquaresProblems,
270
276
check out the [ NonlinearSolve.jl tutorials!] ( https://docs.sciml.ai/NonlinearSolve/stable/tutorials/getting_started/ ) .
271
277
@@ -293,8 +299,8 @@ to see the problem is not equation 2 but other equations in the system. Meanwhil
293
299
some of the conditions:
294
300
295
301
``` @example init
296
- iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
297
- [x => 1, y => 0.0, D(y) => 0.0, λ => 0], [g => 1], guesses = [λ => 1])
302
+ iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
303
+ [x => 1, y => 0.0, D(y) => 0.0, λ => 0], [g => 1], guesses = [λ => 1])
298
304
```
299
305
300
306
gives a NonlinearLeastSquaresProblem which can be solved:
@@ -309,10 +315,9 @@ sol.resid
309
315
310
316
In comparison, if we have a well-conditioned system:
311
317
312
-
313
318
``` @example init
314
- iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
315
- [x => 1, y => 0.0], [g => 1], guesses = [λ => 1])
319
+ iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
320
+ [x => 1, y => 0.0], [g => 1], guesses = [λ => 1])
316
321
```
317
322
318
323
notice that we instead obtained a NonlinearSystem. In this case we have to use
@@ -348,7 +353,7 @@ by initializing the derivatives to zero:
348
353
349
354
``` @example init
350
355
prob = ODEProblem(simpsys, [D(x) => 0.0, D(y) => 0.0], tspan, guesses = [x => 1, y => 1])
351
- sol = solve(prob, Tsit5(), abstol= 1e-16)
356
+ sol = solve(prob, Tsit5(), abstol = 1e-16)
352
357
```
353
358
354
359
Notice that this is a "numerical zero", not an exact zero, and thus the solution will leave the
@@ -370,4 +375,4 @@ sol[α * x - β * x * y]
370
375
371
376
``` @example init
372
377
plot(sol)
373
- ```
378
+ ```
0 commit comments