diff --git a/semester3/numcs/numcs-summary.pdf b/semester3/numcs/numcs-summary.pdf index 82422ad..8ef6128 100644 Binary files a/semester3/numcs/numcs-summary.pdf and b/semester3/numcs/numcs-summary.pdf differ diff --git a/semester3/numcs/parts/05_curve-fitting/02_non-linear/02_gauss-newton.tex b/semester3/numcs/parts/05_curve-fitting/02_non-linear/02_gauss-newton.tex index a42fc5b..ef6c7ce 100644 --- a/semester3/numcs/parts/05_curve-fitting/02_non-linear/02_gauss-newton.tex +++ b/semester3/numcs/parts/05_curve-fitting/02_non-linear/02_gauss-newton.tex @@ -20,60 +20,72 @@ Die Iterationsvorschrift ist gegeben durch: \end{align*} \begin{code}{python} -import numpy as np +def newton(x: np.ndarray, GF, HF, tol=1e1-6, maxIter=50): + """ Newton method requires the Gradient & Hessian, more accurate """ + s = np.linalg.solve(HF(x), GF(x)) + x -= s + k = 1 + while np.linalg.norm(s) > tol * np.linalg.norm(x) and k < maxIter: + s = np.linalg.solve(HF(x), GF(x)) + x -= s + k += 1 + return x, k +\end{code} -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 +\begin{code}{python} +def gauss_newton(x: np.ndarray, F, DF, tol=1e-6, maxIter=50): + """ Gauss-Newton algorithm to solve non-linear problem, needs 'only' the Jacobian """ + s = np.linalg.lstsq(DF(x), F(x))[0] + x = x-s + k = 1 + while np.linalg.norm(s) > tol * np.linalg.norm(x) and k < maxIter: + s = np.linalg.lstsq(DF(x), F(x))[0] + x = x-s + k += 1 + return x, k \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). +\inlineex Wir haben zwei Modellfunktionen, $F_1(t) = a_1 + b_1 e^{-c_1 t}$ und $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 - +# Data given 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 +heat = np.array([24.34, 18.93, 17.09, 16.27, 15.97, 15.91]) +cool = np.array([9.66, 18.8, 22.36, 24.07, 24.59, 24.91]) + +# Functions & Jacobian +F_1 = lambda a: a[0] + a[1]*np.exp( -a[2]*t ) - heat +F_2 = lambda a: a[0] - a[1]*np.exp( -a[2]*t ) - cool -# define the corresponding Jacobi matrices def J_1(a): - mat = np.zeros((n, 3)) + J = 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 + J[k, 0] = 1 + J[k, 1] = np.exp( -t[k] * a[2] ) + J[k, 2] = -t[k]*a[1] * np.exp( -t[k]*a[2] ) + return J def J_2(a): - mat = np.zeros((n, 3)) + J = 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 + J[k, 0] = 1 + J[k, 1] = -np.exp( -t[k]*a[2] ) + J[k, 2] = t[k]*a[1] * np.exp( -t[k]*a[2] ) + return J -# guess starting vector -x_1 = np.array([10.0, 5.0, 0.0]) -x_2 = np.array([30.0, 10.0, 0.0]) +# Start vectors (Guessed) +x_1 = np.array([10, 5, 0]) +x_2 = np.array([30, 10, 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) +# Coefficients via Gauss-Newton +a_1, it_1 = gauss_newton(x_1, F_1, J_1) +a_2, it_2 = gauss_newton(x_2, F_2, J_2) + +# The final functions found via Gauss-Newton +f_1 = lambda x: a_1[0] + a_1[1]*np.exp( -a_1[2]*x ) +f_2 = lambda x: a_2[0] - a_2[1]*np.exp( -a_2[2]*x ) \end{code}