diff --git a/week5/ex1.py b/week5/ex1.py index e50b779..cd07c21 100755 --- a/week5/ex1.py +++ b/week5/ex1.py @@ -42,32 +42,38 @@ def mk_phi_rk4(f): return phi_rk4 - # Multi Step #----------- -def mk_phi_AB3(f, phi_short): +def mk_phi_AB3(f, phi_short = None): steps = 3 def phi_AB3(x, y, h): - if len(x) <= steps: - return phi_short(x, y, h) + if len(x) < steps: + if phi_short is None: + raise ValueError("This function needs at least {} steps, {} given.".format(steps, len(x))) + else: + return phi_short(x, y, h) else: return ( 23*f(x[-1], y[-1]) - 16*f(x[-2], y[-2]) + 5*f(x[-3], y[-3]) )/12 return phi_AB3 -def mk_phi_AB4(f, phi_short): +def mk_phi_AB4(f, phi_short = None): steps = 4 def phi_AB4(x, y, h): - if len(x) <= steps: - return phi_short(x, y, h) + if len(x) < steps: + if phi_short is None: + raise ValueError("This function needs at least {} steps, {} given.".format(steps, len(x))) + else: + return phi_short(x, y, h) else: return ( 55*f(x[-1], y[-1]) - 59*f(x[-2], y[-2]) + 37*f(x[-3], y[-3]) -9*f(x[-4],y[-4]))/24 return phi_AB4 + # Integrator # ############## -def integrator(x, y_0, phi): +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) @@ -80,8 +86,14 @@ def integrator(x, y_0, phi): y = np.zeros((N,M), dtype=np.float64) y[0] = y_0 + j = 0 + if y_i is not None: + for i, _ in enumerate(y_i): + y[i+1] = y_i[i] + j = i + h = x[1]-x[0] - for i in range(1,N): + for i in range(j+1, N): y[i] = y[i-1] + h*phi(x[:i], y[:i], h) return y @@ -104,7 +116,6 @@ def test_integrator_singlestep(): exact_sol = ODEF_sol N = 2e1 - M = 1 # Domain on which to do the integration x = np.linspace(0, 1, N) @@ -140,26 +151,35 @@ def test_integrator_multistep(): exact_sol = ODEF_sol N = 2e1 - M = 1 # Domain on which to do the integration x = np.linspace(0, 1, N) # Calculate Initial Values using RK4 y_0 = 0.5 - schemes = [ - ["AB3", mk_phi_AB3], - ["AB4", mk_phi_AB4], - ] + y_init = integrator(x[:4], y_0, mk_phi_rk4(test_func)) + + #Name, func, steps + multi_schemes = [ + ["AB3", mk_phi_AB3, 3], + ["AB4", mk_phi_AB4, 4], + ] # Show Plot from matplotlib import pyplot pyplot.subplots() - for name, func in schemes: - pyplot.plot(x, integrator(x, y_0, func(test_func, mk_phi_rk4(test_func))), '--o', label=name) - + # Exact Solution pyplot.plot(x, exact_sol(x), '-', label="Exact Solution") + + # Plot the schemes + for name, func, steps in multi_schemes: + try: + pyplot.plot(x, integrator(x, y_0, func(test_func), y_init[:steps]), '--o', label=name) + except Exception as e: + print(e) + pass + pyplot.xlabel("x") pyplot.ylabel("y") @@ -167,7 +187,65 @@ def test_integrator_multistep(): pyplot.show() +def plot_integration_errors(): + test_func = ODEF + exact_sol = ODEF_sol + + # Domain on which to do the integration + x = np.arange(0, 2.001, step=0.02) + + # Calculate Initial Values using RK4 + y_0 = 0.5 + y_init = integrator(x[:4], y_0, mk_phi_rk4(test_func)) + + single_schemes = [ + ["Euler", mk_phi_euler], + ["Collatz", mk_phi_euler_mod], + ["Heun", mk_phi_heun], + ["RK4", mk_phi_rk4], + + ] + multi_schemes = [ + ["AB3", mk_phi_AB3, 3], + ["AB4", mk_phi_AB4, 4], + ] + + # Show Plot + from matplotlib import pyplot + pyplot.subplots() + + # Pre calculate the exact solution + exact = exact_sol(x) + exact = np.reshape(exact, (-1,1)) + + # Plot Single Step Schemes + for name, func in single_schemes: + try: + pyplot.plot(x, np.abs(exact - integrator(x, y_0, func(test_func))), '--o', label=name) + except Exception as e: + print(e) + pass + + # Plot Multi Step Schemes + for name, func, steps in multi_schemes: + try: + pyplot.plot(x, np.abs(exact - integrator(x, y_0, func(test_func), y_init[:steps])), '--o', label=name) + except Exception as e: + print(e) + pass + + pyplot.xlabel("x") + pyplot.ylabel("absolute error $|\\bar{y} - y|$") + + pyplot.legend() + + pyplot.show() + + if __name__ == "__main__": np.random.seed(0) #test_integrator_singlestep() - test_integrator_multistep() + #test_integrator_multistep() + + plot_integration_errors() +