ex7: problem 10

exercise7/main.typ

@@ -45,3 +45,13 @@ this document was created using
 = problem 8
 
 #include "problem8.typ"
+
+#pagebreak()
+
+= problem 10
+
+#include "problem10.typ"
+
+#line(length: 100%)
+
+I think I understand at least some numerics and ODEs now... 😅

exercise7/problem10.typ

@@ -0,0 +1,116 @@
+#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.
+