117 lines
2.9 KiB
Typst
117 lines
2.9 KiB
Typst
#import "@preview/physica:0.9.6": *
|
|
|
|
consider
|
|
$
|
|
cases(
|
|
y' = -2 t y^2,
|
|
y(0) = 1
|
|
)
|
|
$
|
|
|
|
== a)
|
|
|
|
we can solve for the exact solution by integrating both sides
|
|
$
|
|
dv(y, t) & = -2 t y^2 \
|
|
1/y^2 dd(y) & = -2 t dd(t) \
|
|
integral y^(-2) dd(y) & = integral -2 t dd(t) \
|
|
-1 / y & = -t^2 + C \
|
|
y & = 1 / (t^2 + C) \
|
|
$
|
|
and since $y(0) = 1$ we obtain
|
|
$
|
|
y(0) = 1/(0^2 + C) = 1/C = 1 => C = 1
|
|
$
|
|
thus
|
|
$
|
|
y(t) = 1/(t^2 + 1).
|
|
$
|
|
|
|
plugging in $t = 0.4$ yields
|
|
$
|
|
y(0.4) = y(2/5) & = 1/((2/5)^2 + 1) \
|
|
& = 1 / (4/25 + 1) \
|
|
& = 25 / (4 + 25) \
|
|
& = 25 / 29 \
|
|
& = 0.86206896551724
|
|
$
|
|
as expected.
|
|
|
|
|
|
== b)
|
|
|
|
recall that euler's method is $y_(n+1) = y_n + h f(t_n, y_n)$, and let $f(t, y)
|
|
= -2 t y^2$.
|
|
|
|
then, with a step size of $h = 0.1$ we get
|
|
#let euler = 0.941584 + 0.1 * (-2 * 0.3 * 0.941584 * 0.941584)
|
|
$
|
|
y(0) & = 1, \
|
|
y(0.1) & = 1 + 0.1 dot (-2 dot 0 dot 1^2) = 1, \
|
|
y(0.2) & = 1 + 0.1 dot (-2 dot 0.1 dot 1^2) = 1 - 0.02 = 0.98, \
|
|
y(0.3) & = 0.98 + 0.1 dot (-2 dot 0.2 dot 0.98^2) = 0.941584, \
|
|
y(0.4) & = #calc.round(euler, digits: 6)
|
|
$
|
|
after four steps of euler's method.
|
|
|
|
the error is
|
|
$
|
|
e_4 := abs(y_4 - y(0.4)) = #calc.round(calc.abs(euler - 25 / 29), digits: 6).
|
|
$
|
|
|
|
|
|
== c)
|
|
|
|
recall heun's method is
|
|
$
|
|
tilde(y)_(n+1) & = y_n + h f(t_n, y_n) quad #text(gray)[euler's method] \
|
|
y_(n+1) & = y_n + h/2 [f(t_n, y_n) + f(t_(n+1), tilde(y)_(n+1))]
|
|
$
|
|
|
|
thus with $h = 0.2$, two steps of heun's method looks like
|
|
#let heun = 0.911592 + 0.1 * (-0.39204 - 2 * 0.4 * 0.911592 * 0.911592)
|
|
$
|
|
y(0) & = 1 \
|
|
tilde(y)(0.2) & = 1 + 0.2 dot (-2 dot 0 dot 1^2) = 1 \
|
|
y(0.2) & = 1 + 0.1 [0 - 2 dot 0.2 dot 1^2] = 0.99 \
|
|
tilde(y)(0.4) & = 0.99 + 0.2 dot (-2 dot 0.2 dot 0.99^2) = 0.911592 \
|
|
y(0.4) & = 0.911592 + 0.1 dot (-0.39204 - 2 dot 0.4 dot 0.911592^2) \
|
|
& = #calc.round(heun, digits: 6)
|
|
$
|
|
|
|
this yields an error of
|
|
$
|
|
e_2 := #calc.round(25 / 29 - heun, digits: 6)
|
|
$
|
|
which is greater than for euler. curious.
|
|
|
|
|
|
== d)
|
|
|
|
$
|
|
k_1 & = f(t_n, y_n) = f(0, 1) = 0 \
|
|
k_2 & = f(t_n + h/2, y_n + h/2 k_1) = -2 dot 0.2 dot 1^2 = -0.1 \
|
|
k_3 & = f(t_n + h/2, y_n + h/2 k_2) = -2 dot 0.2 dot 0.98^2 = -0.38416 \
|
|
k_4 & = f(t_n + h, y_n + h k_3) = -2 dot 0.4 dot (1 - 0.4 dot 0.38416) = -0.6770688
|
|
$
|
|
thus we obtain
|
|
#let rk4 = 1 + 0.4 / 6 * (0 + 2 * (-0.1) + 2 * (-0.38416) - 0.6770688)
|
|
$
|
|
y_(n+1) & = y_n + h/6 (k_1 + 2 k_2 + 2 k_3 + k_4) \
|
|
& = 1 + 0.4/6 (0 + 2 (-0.1) + 2 (-0.38416) - 0.6770688 ) \
|
|
& = #calc.round(rk4, digits: 6)
|
|
$
|
|
|
|
which yields the error
|
|
$
|
|
e_1 := #{ calc.round(calc.abs(25 / 29 - rk4), digits: 6) }
|
|
$
|
|
|
|
which is also worse than euler, suprpsingly, but impressive for only a single
|
|
step.
|
|
|
|
we can conclude that euler performed best, but it was given a much finer
|
|
granularity, requiring four explicit steps to get to something slightly better
|
|
than that which rk4 achieved in a single step.
|
|
|