This nimib example document shows how to:
- apply (or not) a custom style (latex.css); for default style click here
- using latex rendering of equations (thanks to katex)
- having as output an html table (using standard markdown converter to table)
The document itself shows how to use numericalnim
to integrate an ODE.
Using NumericalNim
Example of usage of numericalnim.
import math, numericalnim
ODE
We want to solve the IVP (Initial Value Problem) for this linear first order differential equation:
$$y' - y = - \frac{1}{2} e^{\frac{t}{2}} \sin(5t) + 5e^{\frac{t}{2}}\cos(5t) \qquad y(0) = 0$$
This ODE has an analytical solution:
$$y(t) = e^{\frac{t}{2}}\sin(5t)$$
We want to find the approximation to the solution and compare it with the analytical solution at $t=5$.
We will use two fixed timestep methods to find the solution:
- Heun's 2nd order method (
Heun2)
- classic 4th order Runge-Kutta (
RK4)
As timestamps we will use $h=0.1, 0.05, 0.01$.
First we translate the ODE as $y'=f(y,t)$ with $f$:
proc f(t, y: float, ctx: NumContext[float, float]): float =
y - 0.5*exp(0.5*t)*sin(5*t) + 5*exp(0.5*t)*cos(5*t)
proc y(t: float): float =
exp(0.5*t)*sin(5*t)
let y0 = 0.0
var ctx = newNumContext[float, float]()
echo "y'(0): ", f(0, y0, ctx)
let y5 = y(5)
echo "y(5): ", y5
y'(0): 5.0
y(5): -1.612374396254655
We want the solution to be available for every point in $[0, 5]$ with timestep $h=0.05$
let tspan = arange(0.0, 5.0, 0.05, includeEnd=true)
echo tspan[0 .. 2], " . . . ", tspan[^3 .. ^1]
@[0.0, 0.05, 0.1] . . . @[4.9, 4.95, 5.0]
We compute the solution according to our selected 2 methods and 3 timesteps:
let
options1 = newODEoptions(dt = 0.1)
options2 = newODEoptions(dt = 0.05)
options3 = newODEoptions(dt = 0.01)
let
(t1hn, y1hn) = solveODE(f, y0, tspan, options = options1, integrator = "Heun2")
(t2hn, y2hn) = solveODE(f, y0, tspan, options = options2, integrator = "Heun2")
(t3hn, y3hn) = solveODE(f, y0, tspan, options = options3, integrator = "Heun2")
(t1rk, y1rk) = solveODE(f, y0, tspan, options = options1, integrator = "RK4")
(t2rk, y2rk) = solveODE(f, y0, tspan, options = options2, integrator = "RK4")
(t3rk, y3rk) = solveODE(f, y0, tspan, options = options3, integrator = "RK4")
assert t1hn == tspan
assert t2hn == tspan
assert t3hn == tspan
assert t1rk == tspan
assert t2rk == tspan
assert t3rk == tspan
echo "Heun2:"
echo fmt"{y1hn[^1]:.7f}"
echo fmt"{y2hn[^1]:.7f}"
echo fmt"{y3hn[^1]:.7f}"
echo "RK4:"
echo fmt"{y1rk[^1]:.7f}"
echo fmt"{y2rk[^1]:.7f}"
echo fmt"{y3rk[^1]:.7f}"
echo "Analytical:"
echo fmt"{y5:.7f}"
Heun2:
-1.1273196
-1.4821550
-1.6068528
RK4:
-1.6131927
-1.6124233
-1.6123745
Analytical:
-1.6123744
As expected Heun is not very accurate even with smaller timesteps, while Runge-Kutta is very reliable even with bigger timesteps.
To compute percentage error of each method at $y(5)$ I will use:
proc pe(yApprox: float): string = fmt"{(100.0*abs(yApprox - y5) / abs(y5)):.3f}%"
echo pe y1hn[^1]
30.083%
The following table is built as a Markdown table in a nbText block.
Table 1. Percentage errors
| Timestep |
Heun2 |
RK4 |
| $h=0.1$ |
30.083% |
0.051% |
| $h=0.05$ |
8.076% |
0.003% |
| $h=0.01$ |
0.342% |
0.000% |
and here is the code for this Markdown block:
fmt"""
The following table is built as a Markdown table in a `nbText` block.
**Table 1.** Percentage errors
Timestep | Heun2 | RK4
---------|--------------:|-------------:
$h=0.1$ | {pe y1hn[^1]} | {pe y1rk[^1]}
$h=0.05$ | {pe y2hn[^1]} | {pe y2rk[^1]}
$h=0.01$ | {pe y3hn[^1]} | {pe y3rk[^1]}
and here is the code for this Markdown block:
"""
import nimib, strformat, strutils
nbInit
nb.useLatex
let
filename_default_style = nb.filename.replace(".html", "_default_style.html")
default_style = nb.context["stylesheet"]
nb.context["stylesheet"] = """<link rel="stylesheet" href="https://latex.now.sh/style.css">"""
proc introText(defaultStyle = false): string =
let otherStyle = if defaultStyle:
fmt"; _you are looking at default style, for custom style [click here]({(nb.filename.AbsoluteFile).relPath})_"
else:
fmt"; _for default style [click here]({(filename_default_style.AbsoluteFile).relPath})_"
fmt"""
> This nimib example document shows how to:
> - apply (or not) a custom style ([latex.css](https://latex.now.sh/)){otherStyle}
> - using latex rendering of equations (thanks to [katex](https://katex.org/))
> - having as output an html table (using standard markdown converter to table)
>
> The document itself shows how to use [numericalnim](https://github.com/hugogranstrom/numericalnim)
> to integrate an ODE.
"""
nbText: introText()
var blockIntroText = nb.blk
nbText: fmt"""# Using NumericalNim
Example of usage of [numericalnim](https://github.com/hugogranstrom/numericalnim).
"""
nbCode:
import math, numericalnim
nbText: """
## ODE
### Example 3 from [Paul's Online Notes](https://tutorial.math.lamar.edu/classes/de/eulersmethod.aspx)
We want to solve the IVP (Initial Value Problem) for this linear first order differential equation:
$$y' - y = - \frac{1}{2} e^{\frac{t}{2}} \sin(5t) + 5e^{\frac{t}{2}}\cos(5t) \qquad y(0) = 0$$
This ODE has an analytical solution:
$$y(t) = e^{\frac{t}{2}}\sin(5t)$$
We want to find the approximation to the solution and compare it with the analytical solution at $t=5$.
We will use two fixed timestep methods to find the solution:
- Heun's 2nd order method (`Heun2`)
- classic 4th order Runge-Kutta (`RK4`)
As timestamps we will use $h=0.1, 0.05, 0.01$.
"""
nbText: "First we translate the ODE as $y'=f(y,t)$ with $f$:"
nbCode:
proc f(t, y: float, ctx: NumContext[float, float]): float =
y - 0.5*exp(0.5*t)*sin(5*t) + 5*exp(0.5*t)*cos(5*t)
proc y(t: float): float =
exp(0.5*t)*sin(5*t)
let y0 = 0.0
var ctx = newNumContext[float, float]()
echo "y'(0): ", f(0, y0, ctx)
let y5 = y(5)
echo "y(5): ", y5
nbText: "We want the solution to be available for every point in $[0, 5]$ with timestep $h=0.05$"
nbCode:
let tspan = arange(0.0, 5.0, 0.05, includeEnd=true)
echo tspan[0 .. 2], " . . . ", tspan[^3 .. ^1]
nbText: "We compute the solution according to our selected 2 methods and 3 timesteps:"
nbCode:
let
options1 = newODEoptions(dt = 0.1)
options2 = newODEoptions(dt = 0.05)
options3 = newODEoptions(dt = 0.01)
let
(t1hn, y1hn) = solveODE(f, y0, tspan, options = options1, integrator = "Heun2")
(t2hn, y2hn) = solveODE(f, y0, tspan, options = options2, integrator = "Heun2")
(t3hn, y3hn) = solveODE(f, y0, tspan, options = options3, integrator = "Heun2")
(t1rk, y1rk) = solveODE(f, y0, tspan, options = options1, integrator = "RK4")
(t2rk, y2rk) = solveODE(f, y0, tspan, options = options2, integrator = "RK4")
(t3rk, y3rk) = solveODE(f, y0, tspan, options = options3, integrator = "RK4")
assert t1hn == tspan
assert t2hn == tspan
assert t3hn == tspan
assert t1rk == tspan
assert t2rk == tspan
assert t3rk == tspan
echo "Heun2:"
echo fmt"{y1hn[^1]:.7f}"
echo fmt"{y2hn[^1]:.7f}"
echo fmt"{y3hn[^1]:.7f}"
echo "RK4:"
echo fmt"{y1rk[^1]:.7f}"
echo fmt"{y2rk[^1]:.7f}"
echo fmt"{y3rk[^1]:.7f}"
echo "Analytical:"
echo fmt"{y5:.7f}"
nbText: "As expected Heun is not very accurate even with smaller timesteps, while Runge-Kutta is very reliable even with bigger timesteps."
nbText: "To compute percentage error of each method at $y(5)$ I will use:"
nbCode:
proc pe(yApprox: float): string = fmt"{(100.0*abs(yApprox - y5) / abs(y5)):.3f}%"
echo pe y1hn[^1]
nbTextWithCode: fmt"""
The following table is built as a Markdown table in a `nbText` block.
**Table 1.** Percentage errors
Timestep | Heun2 | RK4
---------|--------------:|-------------:
$h=0.1$ | {pe y1hn[^1]} | {pe y1rk[^1]}
$h=0.05$ | {pe y2hn[^1]} | {pe y2rk[^1]}
$h=0.01$ | {pe y3hn[^1]} | {pe y3rk[^1]}
and here is the code for this Markdown block:
"""
let mdCode = nb.blk.code
nbCode:
discard
nb.blk.code = mdCode
nbSave
blockIntroText.output = introText(default_style=true)
nb.filename = filename_default_style
nb.context["stylesheet"] = default_style
nbSave