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.

#### 151 lines 4.1 KiB Raw Permalink Blame History

 ```#!/usr/bin/env python3 ``` ``` ``` ```import numpy as np ``` ```from matplotlib import pyplot ``` ``` ``` ```from solvers import * ``` ``` ``` ```def linear_equation_system_1(): ``` ``` A = [[ 10, -1, 2, 0 ], ``` ``` [ -1, 11, -1, 3 ], ``` ``` [ 2, -1, 10, -1 ], ``` ``` [ 0, 3, -1, 8 ]] ``` ``` ``` ``` b = [ 6, 25, -11, 15 ] ``` ``` ``` ``` return A, b ``` ``` ``` ```def linear_equation_system_2(): ``` ``` A = [[ 0.2, 0.1, 1.0, 1.0, 0.0], ``` ``` [ 0.1, 4.0, -1.0, 1.0, -1.0], ``` ``` [ 1.0, -1.0, 60.0, 0.0, -2.0], ``` ``` [ 1.0, 1.0, 0.0, 8.0, 4.0], ``` ``` [ 0.0, -1.0, -2.0, 4.0, 700.0]] ``` ``` ``` ``` b = [ 1, 2, 3, 4, 5 ] ``` ``` ``` ``` return A, b ``` ``` ``` ``` ``` ```def get_solver_stats(A = None, b = None, epsilons = None, solvers = None): ``` ``` import timeit ``` ``` ``` ``` if A is None: ``` ``` A, b = linear_equation_system_1() ``` ``` ``` ``` if epsilons is None: ``` ``` epsilons = np.array([ 0.1, 0.01, 0.001, 0.0001 ]) ``` ``` ``` ``` if solvers is None: ``` ``` solvers = [ jacobi, steepest_descent, conjugate_gradient ] ``` ``` ``` ``` x_exact = np.linalg.solve(A, b) ``` ``` ``` ``` N = len(solvers) ``` ``` M = len(epsilons) ``` ``` ``` ``` steps = np.zeros((N,M), dtype=np.int) ``` ``` diffs = np.zeros((N,M), dtype=np.float64) ``` ``` times = np.zeros((N,M), dtype=np.float64) ``` ``` ``` ``` for i, func in enumerate(solvers): ``` ``` for j, eps in enumerate(epsilons): ``` ``` start = timeit.timeit() ``` ``` x, k = func(A, b, eps) ``` ``` ``` ``` diffs[i][j] = diff(x_exact, x) ``` ``` steps[i][j] = k ``` ``` times[i][j] = start - timeit.timeit() ``` ``` ``` ``` return steps, diffs, times ``` ``` ``` ```def plot_steps(epsilons, steps, names, fmt = None, show = True): ``` ``` pyplot.subplots() ``` ``` ``` ``` pyplot.title("Performance for different Iterative Methods") ``` ``` pyplot.xlabel("Accuracy \$\epsilon\$") ``` ``` pyplot.ylabel("Steps \$N\$") ``` ``` pyplot.yscale('log') ``` ``` pyplot.xscale('log') ``` ``` ``` ``` if fmt is None: ``` ``` fmt = '-d' ``` ``` ``` ``` for i, name in enumerate(names): ``` ``` if isinstance(fmt, list): ``` ``` pyplot.plot(epsilons, steps[i], fmt[i], label=name) ``` ``` else: ``` ``` pyplot.plot(epsilons, steps[i], fmt, label=name) ``` ``` ``` ``` pyplot.legend() ``` ``` ``` ``` if show: ``` ``` pyplot.show() ``` ``` ``` ```# Main functions ``` ```def main1(): ``` ``` solvers = [ jacobi, steepest_descent, conjugate_gradient ] ``` ``` ``` ``` ## System 1 ``` ``` print("Linear Equation System 1") ``` ``` print(24*"-") ``` ``` A, b = linear_equation_system_1() ``` ``` epsilons = np.logspace(-1, -4, 4) ``` ``` steps, diffs, times = get_solver_stats(A, b, epsilons, solvers) ``` ``` for i, func in enumerate(solvers): ``` ``` print("Name: ", func.__qualname__) ``` ``` print("Eps: ", epsilons) ``` ``` print("Diffs:", diffs[i]) ``` ``` print("Steps:", steps[i]) ``` ``` print("Times:", times[i]) ``` ``` print() ``` ``` ``` ``` plot_steps(epsilons, steps, map(lambda x: x.__name__, solvers)) ``` ``` ``` ``` ``` ```def main2(): ``` ``` solvers = [ jacobi, steepest_descent, conjugate_gradient ] ``` ``` # System 2 ``` ``` print("Linear Equation System 2") ``` ``` print("------------------------") ``` ``` ``` ``` A, b = linear_equation_system_2() ``` ``` epsilons = np.logspace(-1, -8, 8) ``` ``` ``` ``` steps, diffs, times = get_solver_stats(A, b, epsilons, solvers) ``` ``` ``` ``` print("Condition Number:", np.linalg.cond(A)) ``` ``` print(""" ``` ``` Both Steepest Descent and Conjugate Gradient are fast for small accuracies. ``` ``` However, they take many more steps than the Jacobi method when we require larger accuracies. ``` ``` This can be explained by a large Condition Number. ``` ``` """) ``` ``` ``` ``` ``` ``` # Preconditioning of System 2 ``` ``` print("\nConditioned Linear Equation System 2") ``` ``` print( "------------------------------------") ``` ``` D = np.diag(np.array(A, dtype=np.float64)) ``` ``` C = np.diagflat(np.sqrt(np.reciprocal(D))) ``` ``` ``` ``` A_tilde = np.dot(np.dot(C,A), C) ``` ``` b_tilde = np.dot(C,b) ``` ``` ``` ``` cond_steps, cond_diffs, cond_times = get_solver_stats(A_tilde, b_tilde, epsilons, solvers) ``` ``` ``` ``` print("Condition Number:", np.linalg.cond(A_tilde)) ``` ``` print(""" ``` ``` By Preconditioning the Gradient Methods are much faster. ``` ``` This is especially visible for small accuracies in the conjugate gradient method. ``` ``` """) ``` ``` ``` ``` all_steps = np.concatenate((steps, cond_steps), axis=0) ``` ``` fmts = ['b-d', 'g-*', 'y-^', 'b--d', 'g--*', 'y--^'] ``` ``` names = [x.__name__ for x in solvers] + ["Conditioned " + x.__name__ for x in solvers] ``` ``` plot_steps(epsilons, all_steps, names, fmts) ``` ``` ``` ``` ``` ```if __name__ == "__main__": ``` ``` main1() ``` ``` main2() ``` ``` ``` ``` ```