[NumCS] Finish non-linear curve fitting

semester3/numcs/parts/05_curve-fitting/02_non-linear/02_gauss-newton.tex

@@ -18,3 +18,62 @@ Die Iterationsvorschrift ist gegeben durch:
 \begin{align*}
     x^{(k + 1)} = x^{(k)} - s \smallhspace \text{ mit } s := \argmin{z \in \R^n} ||F(x^{(k)}) - DF(x^{(k)})z||^2_2
 \end{align*}
+
+\begin{code}{python}
+import numpy as np
+
+def gauss_newton(start_vec, Func, Jacobian, tolerance):
+    # Start vector has to be chosen intelligently
+    s = np.linalg.lstsq(Jacobian(start_vec), Func(start_vec))[0]
+    start_vec = start_vec - s
+    # now we perform the iteration
+    while np.linalg.norm(s) > tolerance * np.linalg.norm(start_vec):
+        # every time we update x by subtracting s, found with the least square method
+        s = np.linalg.lstsq(Jacobian(start_vec), Func(start_vec))[0]
+        start_vec = start_vec - s
+    return start_vec
+\end{code}
+
+Der Vorteil ist, dass die zweite Ableitung nicht benötigt wird, jedoch ist die Konvergenzordnung niedrieger ($p \leq 2$)
+
+\newpage
+\setLabelNumber{all}{3}
+\inlineex Wir haben zwei Modellfunktionen, $F_1(t) = a_1 + b_1 e^{-c_1 t}$ and $F_2(t) = a_2 - b_2 e^{-c_2 t}$. ($F_1$ ist ein Heizvorgang, $F_2$ ist ein Abkühlvorgang).
+Untenstehender code berechnet die Lösung des nichtlinearen Ausgleichsproblems
+\begin{code}{python}
+import numpy as np
+
+t = np.arange(0, 30, 5); n = len(t)
+curr_heating = np.array([24.34, 18.93, 17.09, 16.27, 15.97, 15.91])
+curr_cooling = np.array([9.66, 18.8, 22.36, 24.07, 24.59, 24.91])
+# define the functions that have to be minimized
+F_1 = lambda a: a[0] + a[1] * np.exp(-a[2] * t) - curr_heating
+F_2 = lambda a: a[0] - a[1] * np.exp(-a[2] * t) - curr_cooling
+
+# define the corresponding Jacobi matrices
+def J_1(a):
+    mat = np.zeros((n, 3))
+    for k in range(n):
+        mat[k, 0] = 1.0
+        mat[k, 1] = np.exp(-t[k] * a[2])
+        mat[k, 2] = -t[k] * a[1] * np.exp(-t[k] * a[2])
+    return mat
+
+def J_2(a):
+    mat = np.zeros((n, 3))
+    for k in range(n):
+        mat[k, 0] = 1.0
+        mat[k, 1] = -np.exp(-t[k] * a[2])
+        mat[k, 2] = t[k] * a[1] * np.exp(-t[k] * a[2])
+    return mat
+
+# guess starting vector
+x_1 = np.array([10.0, 5.0, 0.0])
+x_2 = np.array([30.0, 10.0, 0.0])
+
+# use the Gauss-Newton algorithm declared above
+a_1 = gauss_newton(x_1, F_1, J_1, tolerance=10e-6)
+a_2 = gauss_newton(x_2, F_2, J_2, tolerance=10e-6)
+print("Heating ", a_1)
+print("Cooling ", a_2)
+\end{code}