ericteunis/uni-m.cds-num-met
Archived
1
0
Fork 0
Code Releases Activity

Week6: Wrap up: still the same error

Browse source
This commit is contained in:
Eric Teunis de Boone 2020-03-12 17:27:58 +01:00
parent bda69d945f
commit 160abf9d46
1 changed files with 33 additions and 55 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

74
week6/ex1.py
View file

@ -2,30 +2,23 @@
import numpy as np import numpy as np
def pdeHyperbolic(a, x, t, f, g = None, dtype=np.float64): def pdeHyperbolic(a, x, t, f, g, dtype=np.float64):
""" Solve a Hyperbolic Partial Differential using finite differences. """ """ Solve a Hyperbolic Partial Differential using finite differences. """
m = len(x) # Amount of objects to track m = len(x) # Amount of objects to track
n = len(t) # Length of Time Vector n = len(t) # Length of Time Vector
# Determine stepsizes # Determine stepsizes
h = 1
if m > 1:
h = x[0] - x[1] h = x[0] - x[1]
k = 1
if n > 1:
k = t[0] - t[1] k = t[0] - t[1]
λ_sq = (a*k/h)**2 λ_sq = (a*k/h)**2
# Create array to hold the solution # Create array to hold the solution
w = np.zeros((n,m), dtype=dtype) w = np.zeros((n,m), dtype=dtype)
# Create finite difference matrix # 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) 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 # Initialise first two timesteps
w[0] = f(x, t[0]) w[0] = f(x, t[0])
@ -38,16 +31,13 @@ def pdeHyperbolic(a, x, t, f, g = None, dtype=np.float64):
return w return w
def test_pdeHyperbolic_case1(x_steps=1e2, t_steps=1e2, max_x=1, max_t=1):
def test_pdeHyperbolic_case1(m=1e2, n=1e3,l=1, T=1):
a = 1 # from the Schroedinger Equation a = 1 # from the Schroedinger Equation
# Setup spatial and time grids # Setup spatial and time grids
x = np.linspace(0, l, m) x = np.linspace(0, max_x, x_steps)
t = np.linspace(0, T, n) t = np.linspace(0, max_t, t_steps)
# Boundary conditions # Boundary conditions
def f(x,t): def f(x,t):
@ -60,43 +50,18 @@ def test_pdeHyperbolic_case1(m=1e2, n=1e3,l=1, T=1):
sol = pdeHyperbolic(a, x, t, f, g) sol = pdeHyperbolic(a, x, t, f, g)
# Plot it with the exact solution # 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)) exact_f = 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 plot_animation(x, sol, func=exact_f)
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() def test_pdeHyperbolic_case2(x_steps=2e2, t_steps=4e2, max_x=1, max_t=1):
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 a = 1 # from the Schroedinger Equation
# Setup spatial and time grids # Setup spatial and time grids
x = np.linspace(0, l, m) x = np.linspace(0, max_x, x_steps)
t = np.linspace(0, T, n) t = np.linspace(0, max_t, t_steps)
# Boundary conditions # Boundary conditions
def f(x,t): def f(x,t):
@ -108,6 +73,11 @@ def test_pdeHyperbolic_case2(m=200, n=400, l=1, T=1):
# Solve it # Solve it
sol = pdeHyperbolic(a, x, t, f, g) sol = pdeHyperbolic(a, x, t, f, g)
plot_animation(x, sol)
def plot_animation(x, sol, func=None, interval = 10):
from matplotlib import pyplot from matplotlib import pyplot
from matplotlib import animation as anim from matplotlib import animation as anim
fig, _ = pyplot.subplots() fig, _ = pyplot.subplots()
@ -116,16 +86,24 @@ def test_pdeHyperbolic_case2(m=200, n=400, l=1, T=1):
pyplot.clf() pyplot.clf()
pyplot.grid() pyplot.grid()
pyplot.ylim(-1.5, 1.5) pyplot.ylim(-1.5, 1.5)
pyplot.title('t = {}/{}'.format(i, len(sol)))
#if func is not None:
# pyplot.plot(x, func(x, sol[i]), label='exact')
pyplot.plot(x, sol[i,:], label="iter") pyplot.plot(x, sol[i,:], label="iter")
frames = np.arange(1, n) frames = np.arange(0, len(sol))
myAnim = anim.FuncAnimation(fig, animate, frames, interval = 5 )
try:
myAnim = anim.FuncAnimation(fig, animate, frames, interval = interval )
pyplot.legend() pyplot.legend()
pyplot.show() pyplot.show()
except AttributeError:
# This final error is fugly
pass
if __name__ == "__main__": if __name__ == "__main__":
#test_pdeHyperbolic_case1() test_pdeHyperbolic_case1()
test_pdeHyperbolic_case2() #test_pdeHyperbolic_case2()