#!/usr/bin/env python
# coding: utf-8

# # The Wave Equation
# 

# The wave equation is one of the canonical partial differential equationss with which numericists deal, and an example of a hyperbolic PDE. As you probably know, the equation for a parabola is
# 
# \\[ x^2 = a^2 y^2+b. \\]
#  
# Similar to this, the standard wave equation in one dimension is
# 
# \\[ \frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial ^2 u}{\partial x^2}. \\]
#  
# This is a model for a system radiating information, where the speed of the propagation is \\(c\\), whichand has units of  m/s.
# 
# When \\(c\\) is constant, the equation accepts leftward and rightward travelling wave solutions, \\(f(x+ct)\\) and \\(g(x-ct)\\).

# ## Initial Conditions
# The second order ODE can be rewritten as a pair of coupled first order ODEs,
# \\[ \frac{\partial v}{\partial t} = c^2\frac{\partial ^2 u}{\partial t}, \\]
# \\[ \frac{\partial u}{\partial t} = v. \\]
# where the second equation is a rewritten version of the identity, \\(v=\frac{\partial u}{\partial t}\\). From this, it can be seen that we need to provide initial conditions for both \\(u\\) and \\(v\\), i.e
# \\[ u(x,0)= U(x) \\]
# \\[\frac{\partial u}{\partial t} (x,0)= V(x). \\]
# 
# ## Boundary Conditions
# The boundary conditions must supply enough information to control both leftward and rightward travelling waves, so that we need a boundary condition on every boundary.
# 
# Lets build a finite difference model for a wave equation

# In[44]:


# import relevant modules

import numpy
from scipy import linalg
from bqplot import pyplot
from IPython.display import display
import ipywidgets

# set up our discretization

N = 101
Nt= 100
dx = 10.0/(N-1)
dt=1.0
c=1.0
x = numpy.linspace(-10.0, 10.0, N)

# useful matrices. A handles the new step, B the old.

A = numpy.zeros([2*N, 2*N])
B = numpy.zeros([2*N, 2*N])

for j in range(N):
    # Dirichlet boundary conditions
    A[j,j] = 1.0
    A[N-1,N-1] = 1.0
    A[N+j,N+j] = 1.0
    A[2*N-1,2*N-1] = 1.0
for i in range(1,N-1):
    A[i, i] = 1.0
    A[i, N+i-1] =-dt*c*1.0/dx**2/2.0
    A[i, N+i] = +dt*c*2.0/dx**2/2.0
    A[i, N+i+1] = -dt*c*1.0/dx**2/2.0
    A[N+i, i] = -dt/2.0
    A[N+i, N+i] = 1.0
    B[i, i] = 1.0
    B[i, N+i-1] = dt*c*1.0/dx**2/2.0
    B[i, N+i] = -dt*c*2.0/dx**2/2.0
    B[i, N+i+1] = dt*c*1.0/dx**2/2.0
    B[N+i, i] = dt/2.0
    B[N+i, N+i] = 1.0
    


# initial conditions
    
y = numpy.zeros(2*N)
y[:N]=numpy.exp(-x**2)

Y=[]
    
for _ in range(Nt):
    Y.append(y)
    z=numpy.dot(B,y)
    y = linalg.solve(A, z)

Y.append(y)

pyplot.figure(1)
xsc = pyplot.LinearScale(min=-10,max=10)
ysc = pyplot.LinearScale(min=-1,max=1)
axy = pyplot.Axis(label='u', scale=ysc, orientation='vertical', side='left', grid_lines='solid')
axx = pyplot.Axis(label='x', scale=xsc, grid_lines='solid')
lines = pyplot.Lines(x=x, y=Y[0][:N], scales = {'x': xsc, 'y': ysc})
plt = pyplot.Figure(layout=ipywidgets.Layout(width='auto'), marks=[lines],axes=[axx,axy])
def g(t):
    lines.y = Y[int(t/dt)][:N]
    
w = ipywidgets.interactive(g, t=(0.,1.0*Nt,1));
slider = w.children[0]
output = w.children[-1]

play = ipywidgets.Play(
#     interval=10,
    value=0,
    min=0,
    max=Nt,
    step=1,
    description="Press play",
    disabled=False
)


ipywidgets.jslink((play, 'value'), (slider, 'value'))

display(ipywidgets.HBox([play, slider]))
display(plt)




# In[ ]:





# In[ ]: