Coursework for the Master's course Computation Data Science: Numerical Methods
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
This repo is archived. You can view files and clone it, but cannot push or open issues/pull-requests.

#### 99 lines 2.4 KiB Raw Permalink Blame History

 ```#!/usr/bin/env python3 ``` ``` ``` ```import numpy as np ``` ```from itertools import count as count ``` ``` ``` ```def diff(a, b): ``` ``` return np.amax(np.abs(a-b)) ``` ``` ``` ```def jacobi(A, b, eps, max_iter = None): ``` ``` """ Use the Jacobi Method to solve a Linear System. """ ``` ``` ``` ``` A = np.array(A, dtype=np.float64) ``` ``` b = np.array(b, dtype=np.float64) ``` ``` ``` ``` # Determine Diagonal and Upper and Lower matrices ``` ``` D = np.diag(A) ``` ``` L = -np.tril(A, -1) ``` ``` U = -np.triu(A, 1) ``` ``` ``` ``` D_inv = np.diagflat(np.reciprocal(D)) ``` ``` ``` ``` x_0 = D_inv @ b ``` ``` for i in count(): ``` ``` x_1 = D_inv @ ( L + U) @ x_0 ``` ``` ``` ``` # Are we close enough? ``` ``` if diff(x_0, x_1) < eps: ``` ``` return x_1, i ``` ``` ``` ``` # Running out of iterations ``` ``` if max_iter is not None and max_iter >= i: ``` ``` raise RuntimeError("Did not converge in {} steps".format(max_iter)) ``` ``` ``` ``` # Set values for next loop ``` ``` x_0 = x_1 ``` ``` ``` ```def steepest_descent(A, b, eps, max_iter = None): ``` ``` """ Use Steepest Descent to solve a Linear System. """ ``` ``` A = np.array(A, dtype=np.float64) ``` ``` b = np.array(b, dtype=np.float64) ``` ``` ``` ``` ``` ``` x_0 = np.zeros(len(A), dtype=np.float64) ``` ``` ``` ``` for i in count(): ``` ``` Ax = A @ x_0 ``` ``` v = b - Ax ``` ``` t = np.dot(v,v) / np.dot(v, A @ v ) ``` ``` ``` ``` x_1 = x_0 + t*v ``` ``` ``` ``` # Are we close enough? ``` ``` if diff(x_0, x_1) < eps: ``` ``` return x_1, i ``` ``` ``` ``` # Running out of iterations ``` ``` if max_iter is not None and max_iter >= i: ``` ``` raise RuntimeError("Did not converge in {} steps".format(max_iter)) ``` ``` ``` ``` # Set values for next loop ``` ``` x_0 = x_1 ``` ``` ``` ```def conjugate_gradient(A, b, eps, max_iter = None): ``` ``` """ Use the Conjugate Gradient Method to solve a Linear System. """ ``` ``` A = np.array(A, dtype=np.float64) ``` ``` b = np.array(b, dtype=np.float64) ``` ``` ``` ``` ``` ``` # Setup vectors ``` ``` x_0 = np.zeros(len(A), dtype=np.float64) ``` ``` r_0 = b - A @ x_0 ``` ``` v = r_0.copy() ``` ``` ``` ``` for i in count(): ``` ``` Ax = A @ x_0 ``` ``` Av = A @ v ``` ``` ``` ``` r_0_square = np.dot(r_0, r_0) ``` ``` ``` ``` t = r_0_square / np.dot(v, Av ) ``` ``` x_1 = x_0 + t*v ``` ``` ``` ``` r_1 = r_0 - t * Av ``` ``` ``` ``` s = np.dot(r_1, r_1) / r_0_square ``` ``` v = r_1 + s*v ``` ``` ``` ``` ``` ``` # Are we close enough? ``` ``` if diff(x_0, x_1) < eps: ``` ``` return x_1, i ``` ``` ``` ``` # Running out of iterations ``` ``` if max_iter is not None and max_iter >= i: ``` ``` raise RuntimeError("Did not converge in {} steps".format(max_iter)) ``` ``` ``` ``` # Set values for next loop ``` ``` x_0 = x_1 ``` ``` r_0 = r_1 ``` ``` ```