#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.