diff --git a/week5/ex1.py b/week5/ex1.py index cd07c21..71ac651 100755 --- a/week5/ex1.py +++ b/week5/ex1.py @@ -75,6 +75,7 @@ def mk_phi_AB4(f, phi_short = None): ############## def integrator(x, y_0, phi, y_i = None ): x = np.asarray(x) + if isinstance(y_0, (list,np.ndarray)): y_0 = np.asarray(y_0) else: @@ -236,16 +237,85 @@ def plot_integration_errors(): pyplot.xlabel("x") pyplot.ylabel("absolute error $|\\bar{y} - y|$") + pyplot.title("absolute error between an approach and the exact solution") pyplot.legend() pyplot.show() +def accuracy_per_stepsize( debug=True ): + stepsizes = np.asarray([1e-1, 1e-2, 1e-3, 1e-5]) + + # define the problem to be solved + test_func = lambda x,y: x + exact_sol = lambda x: np.exp(x) + + # set domain and intial value + x_min = 0 + x_max = 2 + y_0 = 1 + + exact_end = exact_sol(x_max) + + # define the schemes to use + schemes = [ + ["Euler", mk_phi_euler, 1], + ["Collatz", mk_phi_euler_mod, 1], + ["Heun", mk_phi_heun, 1], + ["RK4", mk_phi_rk4, 1], + ["AB3", mk_phi_AB3, 3], + ["AB4", mk_phi_AB4, 4], + ] + + # get max steps needed for the multistep schemes + max_steps = 0 + for _, _, steps in schemes: + if steps > max_steps: + max_steps = steps + + # pre calculate for multistep integrations + y_init = np.zeros((len(stepsizes),1))# Note the dimensionality of y_0 + max_steps += 1 + for i, h in enumerate(stepsizes): + x = x_min + np.linspace(0, h * max_steps, max_steps, True) + # Note the dimensionality of y_0 + # TODO: fix up this ugliness + z = integrator(x[:4], y_0, mk_phi_rk4(test_func)) + y_init[i] = z[1,:] + + # plot schemes + from matplotlib import pyplot + pyplot.subplots() + for name, scheme, steps in schemes: + diffs = np.zeros(len(stepsizes)) + if debug: + print("Calculating {}".format(name)) + + for i, h in enumerate(stepsizes): + N = int(abs(x_max-x_min)/h) + x = np.linspace(x_min, x_max, N, True) + y = integrator(x, y_0, scheme(test_func), y_init[:steps]) + diffs[i] = np.abs(exact_end - [-1]) + + pyplot.plot(stepsizes, diffs, '--o', label=name) + + pyplot.xlabel('Stepsize $h$') + pyplot.ylabel('Absolute Error') + pyplot.title('Absolute error at the end of integration') + + pyplot.legend() + + pyplot.show() + + +def pendulum_problem(): + pass + if __name__ == "__main__": np.random.seed(0) - #test_integrator_singlestep() + test_integrator_singlestep() #test_integrator_multistep() - plot_integration_errors() - + #plot_integration_errors() + #accuracy_per_stepsize()