In [ ]:
%pylab inline
In [ ]:
x = array([1., 2.])
print "x = ", x
x
In [ ]:
A = array([[0.,3.],[1.,0.]])
A
In [ ]:
# matrix-vector product:
dot(A,x)
In [ ]:
x = linspace(0,4,5) # creates numpy array
x
In [ ]:
p = x
u = 2 + x**2 # component wise squaring
q = vstack([p,u]) # stack the two vectors vertically
q
In [ ]:
dot(A,q) # matrix-matrix multiply
In [ ]:
# transpose operator:
q.T
In [ ]:
# Note: unlike Matlab, there is no distinction between row and column vectors!
x.T
In [ ]:
# This makes adding columns a bit awkward...
# add rows to the transposed matrix and then transpose back.
# (e.g. for ghost cells want to add a cell to each end with 2 components)
qghost = array([0.,0.])
vstack([qghost, q.T, qghost]).T
In [ ]:
# To solve a linear system:
from numpy.linalg import solve
A = array([[0.,3.],[1.,0.]])
b = array([3,4])
x = solve(A,b)
x
In [ ]:
# To compute eigenvalues and eigenvectors:
from numpy.linalg import eig
s, R = eig(A)
print "Eigenvalues: "
print s
print "Eigenvector matrix: "
print R
print "Note that the eigenvectors might not be sorted and the vectors are scaled to have unit 2-norm"
This function approximates the solution from time $t_0 = 0$ to some final time tfinal by taking nsteps time steps with the Lax-Wendroff method, applied to a general linear equation $q_t + Aq_x = 0$. The matrix $A$ (as a numpy array) is one input argument.
On input, x is an array of cell centers and q0 should be an array with q0.shape = (meqn,mx) where meqn is the number of equations in the system (2 for acoustics) and mx is the number of grid cells. q0 should contain the initial cell averages at time 0.
In [ ]:
def lax_wendroff(A,q0,x,tfinal,nsteps):
from numpy.linalg import eig
dt = float(tfinal)/nsteps
dx = x[1] - x[0]
(s,R) = eig(A) # compute eigenvalues and eigenvector matrix
print "Eigenvalues: ",s
cfl = abs(s).max() * dt/dx
print "Courant number is ",cfl
qghost = array([0.,0.])
qn = vstack([qghost, q0.T, qghost]).T # add a ghost cell on each end
mx = len(x) # number of grid cells
for n in range(nsteps):
qn[:,0] = qn[:,mx]
qn[:,mx+1] = qn[:,1]
qnp = zeros(qn.shape) # initialize array
for i in range(1,mx+1):
dq = qn[:,i+1] - qn[:,i-1]
dq2 = qn[:,i+1] - 2*qn[:,i] + qn[:,i-1]
qnp[:,i] = qn[:,i] - 0.5*(dt/dx)*dot(A,dq) + 0.5*(dt/dx)**2 * dot(A,dot(A,dq2))
qn = qnp # for next time step
qfinal = qnp[:,1:(mx+1)] # throw away ghost cells
return qfinal
Set up the computational grid and also a much finer grid for plotting the "exact" solution:
In [ ]:
mx = 100
dx = 1./mx
x = linspace(dx/2, 1.-dx/2, mx) # computational grid
xfine = linspace(0,1,5000) # fine grid for plotting true solution
Note how the periodic boundary conditions are imposed.
Define and plot the initial data:
In [ ]:
def qinit(x):
p = where(abs(x-0.5)<0.1, 3., 2.)
u = zeros(x.shape)
q0 = vstack([p,u])
return q0
q0 = qinit(x)
def plotq(x,q):
subplot(2,1,1)
plot(x,q[0,:], 'bo-', markersize=3)
ylim(1,4)
title('pressure')
subplot(2,1,2)
plot(x,q[1,:], 'bo-', markersize=3)
ylim(-2,2)
title('velocity')
plotq(x,q0)
In [ ]:
# material parameters:
rho = 1.
K = 1.
# coefficient matrix:
A = array([[0.,K], [1/rho,0.]])
print "Sound speed = %g" % sqrt(K/rho)
print "Impedance = %g" % sqrt(K*rho)
Test the Lax-Wendroff method for specific values of the parameters:
In [ ]:
tfinal = 0.2
nsteps = 20
q0 = qinit(x)
q = lax_wendroff(A,q0,x,tfinal,nsteps)
plotq(x,q)
Check that this code works for other parameter values and experiment.
Write a new function godunov that implements the Godunov (upwind) method, using the wave-propagation framework. You can do this for acoustics specifically, but it may be easier to do for general matrices $A$ using the eigen-decomposition that's already computed in the lax_wendroff function. Recall this method takes the form:
$$Q_i^{n+1} = Qi^n - \frac{\Delta t}{\Delta x}\sum{p=1}^m \left((s^p)^+ {\cal W}_{i-1/2}^p
(s^p)^- {\cal W}_{i+1/2}^p \right)$$
where $Q_i^n - Q_{i-1}^n = \sum_{p=1}^m \alpha_{i-1/2}^p r^p = \sum_{p=1}^m {\cal W}_{i-1/2}^p$
Implement a qtrue function to define the true solution and do comparison plots as we did for advection.
In [ ]: