Add fixed-point iteration

book/build.py

@@ -141,7 +141,11 @@ list_of_functions(
             },
             {
                 "title": "Numerical methods",
-                "modules": ["numerics::diff", "numerics::solve"],
+                "modules": [
+                    "numerics::diff",
+                    "numerics::solve",
+                    "numerics::fixed_point",
+                ],
             },
             {
                 "title": "Geometry",

book/src/list-functions-math.md

@@ -395,7 +395,7 @@ fn lcm(a: Scalar, b: Scalar) -> Scalar
 
 ## Numerical methods
 
-Defined in: `numerics::diff`, `numerics::solve`
+Defined in: `numerics::diff`, `numerics::solve`, `numerics::fixed_point`
 
 ### `diff` (Numerical differentiation)
 Compute the numerical derivative of the function \\( f \\) at point \\( x \\) using the central difference method.
@@ -421,6 +421,14 @@ More information [here](https://en.wikipedia.org/wiki/Newton%27s_method).
 fn root_newton<X: Dim, Y: Dim>(f: Fn[(X) -> Y], f_prime: Fn[(X) -> Y / X], x0: X, y_tol: Y) -> X
 ```
 
+### `fixed_point` (Fixed-point iteration)
+Compute the approximate fixed point of a function \\( f: X \rightarrow X \\) starting from \\( x_0 \\), until \\( |f(x) - x| < ε \\).
+More information [here](https://en.wikipedia.org/wiki/Fixed-point_iteration).
+
+```nbt
+fn fixed_point<X: Dim>(f: Fn[(X) -> X], x0: X, ε: X) -> X
+```
+
 ## Geometry
 
 Defined in: `math::geometry`

examples/tests/numerics.nbt

@@ -1,5 +1,6 @@
 use numerics::solve
 use numerics::diff
+use numerics::fixed_point
 
 # Root finding
 
@@ -10,6 +11,12 @@ fn f1_prime(x) = 3 x² - 1
 assert_eq(root_newton(f1, f1_prime, 1, 1e-10), 1.52137970680, 1e-8)
 assert_eq(root_newton(f1, f1_prime, 2, 1e-10), 1.52137970680, 1e-8)
 
+# Fixed point iteration
+let a = 3
+fn f_sqrt3(x: Scalar) = 0.5 * (a / x + x)
+
+assert_eq(fixed_point(f_sqrt3, 1, 1e-10), sqrt(3), 1e-10)
+
 # Differentiation
 
 assert_eq(diff(log, 2.0), 0.5, 1e-5)

numbat/modules/all.nbt

@@ -8,3 +8,4 @@ use extra::algebra
 
 use numerics::diff
 use numerics::solve
+use numerics::fixed_point

numbat/modules/numerics/fixed_point.nbt

@@ -0,0 +1,19 @@
+use core::scalar
+use core::functions
+use core::error
+
+fn _fixed_point<X: Dim>(f: Fn[(X) -> X], x0: X, ε: X, max_iter: Scalar) =
+  if abs(x1 - x0) < ε
+    then x1
+    else
+      if max_iter > 0
+        then _fixed_point(f, x1, ε, max_iter - 1)
+        else error("fixed_point: Exceeded max. number of iterations")
+  where
+    x1 = f(x0)
+
+@name("Fixed-point iteration")
+@url("https://en.wikipedia.org/wiki/Fixed-point_iteration")
+@description("Compute the approximate fixed point of a function $f: X \\rightarrow X$ starting from $x_0$, until $|f(x) - x| < ε$.")
+fn fixed_point<X: Dim>(f: Fn[(X) -> X], x0: X, ε: X) =
+  _fixed_point(f, x0, ε, 100)