Week5: Ex1.1 Almost correct implementation

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()