Mercurial > repos > public > sbplib_julia
comparison src/SimpleTimeSteppers/SimpleTimeSteppers.jl @ 862:a382942e5437 feature/subpackage_simple_timesteppers
Add a few simple timesteppers
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 18 Jan 2022 15:29:54 +0100 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 859:bdc718c38096 | 862:a382942e5437 |
|---|---|
| 1 """ | |
| 2 SimpleTimeSteppers | |
| 3 | |
| 4 Implements a few simple time integration schemes useful for testing the accuracy of spatial discretisations. | |
| 5 """ | |
| 6 module SimpleTimeSteppers | |
| 7 | |
| 8 export rk4 | |
| 9 export rkn4 | |
| 10 export centered_2nd_order_fd | |
| 11 | |
| 12 """ | |
| 13 rk4(F, vₙ, tₙ, Δtₙ) | |
| 14 | |
| 15 Solves ``vₜ = F(v,t)`` using the standard 4th order RK-method. | |
| 16 """ | |
| 17 function rk4(F, vₙ, tₙ, Δtₙ) | |
| 18 k₁ = F(vₙ,tₙ) | |
| 19 k₂ = F(vₙ + Δtₙ/2*k₁, tₙ + Δtₙ/2) | |
| 20 k₃ = F(vₙ + Δtₙ/2*k₂, tₙ + Δtₙ/2) | |
| 21 k₄ = F(vₙ + Δtₙ *k₃, tₙ + Δtₙ ) | |
| 22 | |
| 23 vₙ₊₁ = vₙ + (1/6)*(k₁+2*(k₂+k₃)+k₄)*Δtₙ | |
| 24 | |
| 25 return vₙ₊₁ | |
| 26 end | |
| 27 | |
| 28 """ | |
| 29 rkn4(F, vₙ, v̇ₙ, tₙ, Δtₙ) | |
| 30 | |
| 31 Solves ``vₜₜ = F(v,v̇,t)`` using the Runge-Kutta-Nyström method based on the | |
| 32 standard 4th order RK-method. | |
| 33 """ | |
| 34 function rkn4(F, vₙ, v̇ₙ, tₙ, Δtₙ) | |
| 35 k₁ = F(vₙ, v̇ₙ, tₙ ) | |
| 36 k₂ = F(vₙ + Δtₙ/2*v̇ₙ, v̇ₙ + Δtₙ/2*k₁, tₙ + Δtₙ/2) | |
| 37 k₃ = F(vₙ + Δtₙ/2*v̇ₙ + Δtₙ^2/4*k₁, v̇ₙ + Δtₙ/2*k₂, tₙ + Δtₙ/2) | |
| 38 k₄ = F(vₙ + Δtₙ *v̇ₙ + Δtₙ^2/2*k₂, v̇ₙ + Δtₙ *k₃, tₙ + Δtₙ ) | |
| 39 | |
| 40 vₙ₊₁ = vₙ + Δtₙ*v̇ₙ + Δtₙ^2*(1/6)*(k₁ + k₂ + k₃); | |
| 41 v̇ₙ₊₁ = v̇ₙ + Δtₙ*(k₁ + 2*k₂ + 2*k₃ + k₄)/6; | |
| 42 | |
| 43 return vₙ₊₁, v̇ₙ₊₁ | |
| 44 end | |
| 45 | |
| 46 """ | |
| 47 centered_2nd_order_fd(F, vₙ, vₙ₋₁, Δtₙ) | |
| 48 | |
| 49 Solves ``vₜₜ = F(v,t)`` using the 2nd order discretization | |
| 50 ``(vₙ₊₁ - 2vₙ + vₙ₋₁)/Δt² = F(vₙ,tₙ)``. | |
| 51 """ | |
| 52 function centered_2nd_order_fd(F, vₙ, vₙ₋₁, tₙ, Δt) | |
| 53 vₙ₊₁ = Δt^2*F(vₙ,tₙ) + 2vₙ - vₙ₋₁ | |
| 54 | |
| 55 return vₙ₊₁, vₙ | |
| 56 end | |
| 57 | |
| 58 # TODO: Rewrite without allocations | |
| 59 | |
| 60 | |
| 61 # TODO: Add euler forward and try a local error estimate for testing | |
| 62 end # module |