Skip to content

Commit f39501f

Browse files
committed
Fix energy error plots
Show the change in energy (ΔE), not the energy.
1 parent e2cf401 commit f39501f

File tree

1 file changed

+5
-5
lines changed

1 file changed

+5
-5
lines changed

tutorials/models/01-classical_physics.jmd

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -350,10 +350,10 @@ energy = map(x->E(x...), sol.u)
350350
@show ΔE = energy[1]-energy[end]
351351

352352
#Plot
353-
plot(sol.t, energy, title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")
353+
plot(sol.t, energy .- energy[1], title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")
354354
```
355355

356-
### Symplectic Integration
356+
##### Symplectic Integration
357357

358358
To prevent energy drift, we can instead use a symplectic integrator. We can directly define and solve the `SecondOrderODEProblem`:
359359

@@ -390,17 +390,17 @@ energy = map(x->E(x[3], x[4], x[1], x[2]), sol2.u)
390390
@show ΔE = energy[1]-energy[end]
391391

392392
#Plot
393-
plot(sol2.t, energy, title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")
393+
plot(sol2.t, energy .- energy[1], title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")
394394
```
395395

396-
It's so close to zero it breaks GR! And let's try to use a Runge-Kutta-Nyström solver to solve this. Note that Runge-Kutta-Nyström isn't symplectic.
396+
And let's try to use a Runge-Kutta-Nyström solver to solve this. Note that Runge-Kutta-Nyström isn't symplectic.
397397

398398
```julia
399399
sol3 = solve(prob, DPRKN6());
400400
energy = map(x->E(x[3], x[4], x[1], x[2]), sol3.u)
401401
@show ΔE = energy[1]-energy[end]
402402
gr()
403-
plot(sol3.t, energy, title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")
403+
plot(sol3.t, energy .- energy[1], title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")
404404
```
405405

406406
Note that we are using the `DPRKN6` sovler at `reltol=1e-3` (the default), yet it has a smaller energy variation than `Vern9` at `abs_tol=1e-16, rel_tol=1e-16`. Therefore, using specialized solvers to solve its particular problem is very efficient.

0 commit comments

Comments
 (0)