diff --git a/week6/ex1.py b/week6/ex1.py new file mode 100755 index 0000000..e044ca3 --- /dev/null +++ b/week6/ex1.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python3 + +import numpy as np + +def pdeHyperbolic(a, x, t, f, g = None, dtype=np.float64): + """ Solve a Hyperbolic Partial Differential using finite differences. """ + + m = len(x) # Amount of objects to track + n = len(t) # Length of Time Vector + + # Determine stepsizes + h = 1 + if m > 1: + h = x[0] - x[1] + + k = 1 + if n > 1: + k = t[0] - t[1] + + λ_sq = (a*k/h)**2 + + + # Create array to hold the solution + w = np.zeros((n,m), dtype=dtype) + + # Create finite difference matrix + A = np.diag(m*[2*(1 - λ_sq)], k=0) + np.diag((m-1)*[λ_sq], k=-1) + np.diag((m-1)*[λ_sq], k=1) + print(A) + + # Initialise first two timesteps + w[0] = f(x, t[0]) + w[1] = A@w[0]/2 + k*g(x, t[0]) + + # Calculate for following timesteps + for j in range(2, n-1): + w[j] = A@w[j-1] - w[j-2] + + return w + + + + + +def test_pdeHyperbolic_case1(m=1e2, n=1e3,l=1, T=1): + + a = 1 # from the Schroedinger Equation + + # Setup spatial and time grids + x = np.linspace(0, l, m) + t = np.linspace(0, T, n) + + # Boundary conditions + def f(x,t): + return np.sin(2*np.pi*x) + + def g(x,t): + return 2*np.pi*np.sin(2*np.pi*x) + + # Solve it + sol = pdeHyperbolic(a, x, t, f, g) + + # Plot it with the exact solution + exact_sol = lambda x,t: np.sin(2*np.pi*x)*(np.cos(2*np.pi*t) + np.sin(2*np.pi*t)) + t_high_res = np.linspace(0, T, n * 1e4) + + from matplotlib import pyplot + from matplotlib import animation as anim + fig, _ = pyplot.subplots() + + if False: + for _, x_i in enumerate(x): + pyplot.plot(t_high_res, exact_sol(x_i, t_high_res), label="Exact") + pyplot.plot(t, sol[:,1], label="iter") + + pyplot.grid() + pyplot.xlabel("t") + pyplot.ylabel("w") + + else: + def animate(i): + pyplot.clf() + pyplot.ylim(-2, 2) + pyplot.grid() + pyplot.plot(x, exact_sol(x, t[i]), label="exact") + pyplot.plot(x, sol[i,:], label="iter") + + frames = np.arange(1, n, dtype=np.int) + myAnim = anim.FuncAnimation(fig, animate, frames, interval = 10 ) + + pyplot.legend() + pyplot.show() + +def test_pdeHyperbolic_case2(m=200, n=400, l=1, T=1): + + a = 1 # from the Schroedinger Equation + + # Setup spatial and time grids + x = np.linspace(0, l, m) + t = np.linspace(0, T, n) + + # Boundary conditions + def f(x,t): + return 2*(x < 0.5) -1 + + def g(x,t): + return 0 + + # Solve it + sol = pdeHyperbolic(a, x, t, f, g) + + from matplotlib import pyplot + from matplotlib import animation as anim + fig, _ = pyplot.subplots() + + def animate(i): + pyplot.clf() + pyplot.grid() + pyplot.ylim(-1.5,1.5) + pyplot.plot(x, sol[i,:], label="iter") + + frames = np.arange(1, n) + myAnim = anim.FuncAnimation(fig, animate, frames, interval = 5 ) + + pyplot.legend() + pyplot.show() + + + +if __name__ == "__main__": + #test_pdeHyperbolic_case1() + test_pdeHyperbolic_case2()