#!/usr/bin/env python
'''
Filename: ws9_ex2.py
Author: Denis Jarema
Date: January 10, 2014
Lecture: SciComp
Tutorial: 9
Exercise: 2
-------------------------
'''
import numpy as np
import matplotlib.pyplot as plt
import math
delay = 0.01
# initial setup
timesteps = 100; N = 100; tau = 0.02; h = 1./N; # unstable or stable?
#timesteps = 100; N = 10; tau = 0.02; h = 1./N; # unstable or stable?
#timesteps = 100; N = 100; tau = 0.01; h = 1./N; # unstable or stable?
leftBoundary = 0.; rightBoundary = 0.
points = [ float(k)/N for k in range(0,N+1) ]
# initial conditions
Told = [leftBoundary] +\
[math.exp(-100*(float(k)/N-0.4)**2) for k in range(1,N)] +\
[rightBoundary]
t2bh2 = tau**2/h**2
Tnew = [leftBoundary] +\
[Told[k] + t2bh2/2*(Told[k-1]-2*Told[k]+Told[k+1]) for k in range(1,N)] +\
[rightBoundary]
# time loop
plt.figure()
plt.hold(False)
for t in range(0,timesteps):
Told[1:N] = [t2bh2*(Tnew[k+1]+Tnew[k-1])+2*(1-t2bh2)*Tnew[k]-Told[k] \
for k in range(1,N)]
# swap two lists
Told, Tnew = Tnew, Told
plt.plot(points,Tnew)
plt.axis([0, 1, -1, 1])
plt.pause(delay)
plt.show()