CoCalc Public FilesScalar_Convection_Diffusion.ipynb
Author: Corey Trahan
Views : 114

## Comparison of Traditional MOR with Lagrangian MOR for scalar convection and high-Pe scalar convection-diffusion equations

%%latex

where $f_1 = 1$, $f_2 = 1/Pe$, $w(y,0) = 0.5e^{-(y-0.3)^2/0.05^2}$, $w(0,t) = 0$.

Case 1: $Pe = \infty$ (scalar convection) and Case 2: $Pe = 10^3$ (scalar convection-diffusion)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import sparse
from scipy.sparse import linalg
import random
from scipy import interpolate

In [2]:
%matplotlib inline

In [3]:
# Material Parameters
Pe = 10.0**3

# Case 1 (No Diffusion)
#eps = 0.0
# Case 2 (Convection-Diffusion)
eps = 1.0

f1 = 0.7
f2 = eps*1./Pe

# Homogeneous Dirichlet boundary conditions (compatible with initial condition)
y_l = 0.; y_r = y_l

In [4]:
# spatial interval
xL=0.0; Lx=1.5
# time interval
T=1.0
# initial time
t0=0.
# number of time steps
n_t = 70
# constant time step size
dt =(T-t0)/float(n_t)
DT = dt * np.ones(n_t)

# number of grid points
Nn = 70
# number of cells in grid
Ne = Nn-1
# spatial grid size
dx = (Lx-xL)/float(Ne)
# array of grid points
x = np.linspace(xL,xL+Lx,Nn)

#where the dirichlet boundaries are
# dir_ind = np.array([0,Nn-1])
# out_ind = np.array([],'i')

In [5]:
# Left hand side of 2nd order centrally discretized system
main_diag=np.ones((Nn),)
diagonals = [(1+ 2*dt*f2/dx**2)*main_diag, (0.5*f1*dt/dx - dt*f2/dx**2)*main_diag, -(0.5*f1*dt/dx + f2*dt/dx**2)*main_diag]
offsets = [0, 1, -1]

# resetting coefficient matrix for dirichlet boundary
diagonals[0][0] = 1.0; diagonals[0][-1]= 1.0
diagonals[1][0] = 0.0
diagonals[2][-2]= 0.0

A = scipy.sparse.diags(diagonals, offsets, shape = (Nn,Nn), format = 'csr')

# right hand side for algebraic system
b = np.zeros((Nn,),'d')

In [6]:
# initial data
x0 = x
w0 = 0.5*np.exp(-(x - 0.3)**2/0.05**2)

In [7]:
# Plot initial solution
fig = plt.figure()
ax = fig.gca()
ax.plot(x,w0,'b')
plt.xlabel('x')
plt.ylabel('w')
plt.title('Initial solution')

<matplotlib.text.Text at 0x7fd56e07ced0>

## Generate Snapshots using High Fidelity Model (HFM)

%%latex The discretized High Fidelity Model is given by

In [8]:
# array of exact time point values
tnList = t0*np.ones(n_t+1)
for i in np.arange(1,n_t+1):
tnList[i] = tnList[i-1]+DT[i-1]

t0=tnList[0]; tn = t0; #ii=0;
W = np.zeros((w0.shape[0],len(tnList)),'d')
X = W.copy() # Lagrangian grid calculations
#print X.shape
W[:,0] = w0
X[:,0] = x0;
for nn,t in enumerate(tnList[1:]):
#print "step {0} solving to t={1}".format(nn+1,t)
dt=t-tn #allows for variable time step
X[:,nn+1] = X[:,0] + (nn+1)*dt*f1 # Lagrangian coordinates
b = W[:,nn] #contribution from prev time step
b[0] = y_l # left boundary condition
b[-1] = y_r # right boundary condition
W[:,nn+1] = scipy.sparse.linalg.spsolve(A,b)
tn=t


In [9]:
from ipywidgets import interact

In [10]:
def plot_soln(index_t):
max_w=np.max(W.flat)
min_w=np.min(W.flat)
#print 'Max solution is u={0} at x={1}'.format(max_u,np.where(u == max_u))
# print np.where(u.flat == max_u)
#print 'Min solution is u={0} at x={1}'.format(min_u,np.where(u.flat == min_u))
# print min_u
# print np.where(u.flat == min_u)
fig = plt.figure()
ax = fig.gca()
ax.plot(x,W[:,index_t],'b')
plt.xlabel('x')
plt.ylabel('w')
ax.set_ylim([min_w,max_w])
plt.title('Snapshot at time $t^i=%12.3f$' % (tnList[index_t]))
#leg=ax.legend(['$u^i$','$u^{ex}$'],loc=1)

In [11]:
plot_soln(2*n_t/10)

In [12]:
#interact(plot_soln,index_t=(w.shape[1]-1)/2)


## Perform SVD

Perform SVD for snapshots

$\mathbf{{S}}$ is a diagonal matrix containing the eigenvalues and $\mathbf{{U}}$ is a square matrix containing the eigenvectors of the kernel $\mathbf{{K}}_{t}$,

In [13]:
W_svd  = W.copy()

# The snapshots are centered to have zero mean
#ns_tmp = np.ones(n_t+1,'d')
#w_mean = np.mean(w,axis=1)
#w_svd -= np.outer(w_mean,ns_tmp)

U,S,V  = np.linalg.svd(W_svd,full_matrices=False)
#print np.diag(S)

In [14]:
fig=plt.figure()
ax = fig.gca()
plt.plot(np.log(S),'b*',label='S')
plt.ylabel('singular values of $\mathbf{W}$')

Text(0,0.5,u'singular values of $\\mathbf{W}$')
In [15]:
#truncate solution snapshots
n_k = 3
U_k = U[:,0:n_k]
U_kT= U_k.conj().T
print "Truncating to n_k={0}, \sigma={1}".format(n_k,S[n_k])

#Approximation accuracy estimated from sum of squares of singular values
acc = np.sum(S[0:n_k+1]**2)/np.sum(S**2)
print "Accuracy of reduced representation ={0}".format(acc)

Truncating to n_k=3, \sigma=1.80577869635 Accuracy of reduced representation =0.880742079487
In [16]:
#Rank k approximation
Wpod = U_k.dot(U_kT.dot(W))

In [17]:
def plot_pod(index):
max_w=np.max(W.flat)
min_w=np.min(W.flat)
fig = plt.figure()
ax = fig.gca()
ax.plot(x,W[:,index],'b',x,Wpod[:,index],'r')
plt.xlabel('x')
plt.ylabel('w')
ax.set_ylim([min_w,max_w])
plt.title('solutions at time $t^i=%12.3f$' % (tnList[index]))
leg=ax.legend(['$w^i$','$\hat{w}^i$'],loc=1)

In [18]:
plt.plot(x,W[:,0],'k-.',x,Wpod[:,0],'r-.',x,W[:,n_t/3],'k--',x,Wpod[:,n_t/3],'b--',x,W[:,2*n_t/3],'k-',x,Wpod[:,2*n_t/3],'c-',x,W[:,n_t],'k:',x,Wpod[:,n_t],'m:')
#interact(plot_pod,index=(0,n_t-1))
plt.title('Comparison of truth solution (black) \n and Eulerian MOR (red, blue, cyan and magenta)')

Text(0.5,1,u'Comparison of truth solution (black) \n and Eulerian MOR (red, blue, cyan and magenta)')
In [19]:
def plot_pod_error(index):
diff = W-Wpod
max_err=np.max(diff.flat)
min_err=np.min(diff.flat)
fig = plt.figure()
ax = fig.gca()
ax.plot(x,W[:,index]-Wpod[:,index],'k')
plt.xlabel('x')
plt.ylabel('w')
ax.set_ylim([min_err,max_err])
plt.title('POD error  at time $t^i=%6.3f$' % (tnList[index]))
leg=ax.legend(['$\epsilon_p^i$'],loc=1)

In [20]:
plot_pod_error(n_t/2)
#interact(plot_pod_error,index=(0,n_t-1))


## Perform SVD on Lagrangian coordinate snapshots

Perform SVD for snapshots of Lagrangian coordinates

$\mathbf{\hat{S}}$ is a diagonal matrix containing the eigenvalues and $\mathbf{\hat{U}}$ is a square matrix containing the eigenvectors of the kernel $\mathbf{{K}}_{xt}$,

## Lagrangian MOR

%%latex

The procedure for the case of a linear 1D scalar convection-diffusion equation involves the following steps -

%%latex

Find optimal reduced bases:

Define:

In [21]:
# Perform SVD for Lagrangian coordinates
X_svd = X.copy()

# The snapshots are centered to have zero mean
#ns_tmp = np.ones(n_t+1,'d')
#w_mean = np.mean(w,axis=1)
#w_svd -= np.outer(w_mean,ns_tmp)

Ux,Sx,Vx = np.linalg.svd(X_svd,full_matrices=False)

In [22]:
fig = plt.figure()
ax = fig.gca()
plt.plot(np.log(Sx),'b*',label='Sx')
plt.ylabel('singular values of $\mathbf{X}$')

Text(0,0.5,u'singular values of $\\mathbf{X}$')
In [23]:
# Compute truncated left and right singular matrices for Lagrangian coordinates
nx_k = 3
Ux_k = Ux[:,0:nx_k]
Sx_full = np.diag(Sx)
Vx_k = (Sx_full.dot(Vx))[0:nx_k,:]
Ux_kT = Ux_k.conj().T

#Full Rank approximation of Lagrangian coordinate snapshots
Xpod_full = Ux.dot((Ux.conj().T).dot(X_svd))

#Rank k approximation of Lagrangian coordinate snapshots
Xpod_k = Ux_k.dot(Ux_kT.dot(X_svd))

In [24]:
# Preparing input vector for optimization routine
S_full = np.diag(S)
V_k = (S_full.dot(V))[:n_k,:]
Ux_kflat = Ux_k.flatten('F'); Vx_kflat = Vx_k.flatten('F'); U_kflat = U_k.flatten('F'); V_kflat = V_k.flatten('F')
print 'dim(Ux)={0}, dim(Vx)={1}, dim(Uw)={2}, dim(Vw)={3}'.format(Nn*nx_k, nx_k*(n_t+1), Nn*n_k, n_k*(n_t+1))
print 'dim(Ux)={0}, dim(Vx)={1}, dim(Uw)={2}, dim(Vw)={3}'.format(Ux_k.shape, Vx_k.shape, U_k.shape, V_k.shape)
#print Ux_kflat.shape + Vx_kflat.shape + U_kflat.shape + V_kflat.shape

guess_trunc = np.concatenate((Ux_kflat, Vx_kflat, U_kflat, V_kflat),axis = 0)
guess_full = np.concatenate((U_kflat, V_kflat),axis = 0)
#guess_inp = np.array([Ux_k, Vx_k, U_k, V_k])

print '\nChecking that optimization routine is not indeterminate ....'
print 'Flattened full Lagrangian input guess size = {0}'.format(guess_full.shape)
print 'Flattened truncated Lagrangian input guess size = {0}'.format(guess_trunc.shape)
print 'Flattened residual output size = {0}'.format(len(W.flatten('F')))

dim(Ux)=210, dim(Vx)=213, dim(Uw)=210, dim(Vw)=213 dim(Ux)=(70, 3), dim(Vx)=(3, 71), dim(Uw)=(70, 3), dim(Vw)=(3, 71) Checking that optimization routine is not indeterminate .... Flattened full Lagrangian input guess size = (423,) Flattened truncated Lagrangian input guess size = (846,) Flattened residual output size = 4970

## Function computing the residual error between the true state snapshot basis and the truncated state snapshot basis interpolated on the truncated Lagrangian basis

''' def residual_err2(tilde,full): W_interp1 = W.copy() W_interp2 = W.copy() if full == 1: # FULL Lagrangian basis Uk = np.reshape(tilde[:Nnn_k],(Nn,n_k),'F') Vk = np.reshape(tilde[Nnn_k:],(n_k,n_t+1),'F') ## For full Lagrangian basis X_tilde = Ux.dot(Sx_full.dot(Vx)) #Xpod_k# elif full == 0: # TRUNCATED Lagrangian basis Uxk = np.reshape(tilde[:Nnnx_k],(Nn,nx_k),'F') Vxk = np.reshape(tilde[Nnnx_k:(Nnnx_k + (n_t+1)nx_k)],(nx_k,n_t+1),'F') Uk = np.reshape(tilde[(Nnnx_k + (n_t+1)nx_k):(Nnnx_k + (n_t+1)nx_k + Nnn_k)],(Nn,n_k),'F') Vk = np.reshape(tilde[(Nnnx_k + (n_t+1)nx_k + Nnn_k):],(n_k,n_t+1),'F') ## For truncated Lagrangian basis X_tilde = Uxk.dot(Vxk)

W_tilde = Uk.dot(Vk) #
#W_err = W - W_tilde
for nn,t in enumerate(tnList[1:]):
f1 = interpolate.interp1d(x0, W[:,nn],fill_value="extrapolate")
f2 = interpolate.interp1d(x0, W_tilde[:,nn],fill_value="extrapolate")

W_interp1[:,nn] = f1(X_tilde[:,nn]) # Solution snapshot in Lagrangian coordinates
W_interp2[:,nn] = f2(X_tilde[:,nn]) # Solution snapshot in Lagrangian coordinates
#w_interp[:,nn] = np.interp(x0, X_tilde[:,nn], w_tilde[:,nn])

return (W_interp1 - W_interp2).flatten('F')

''' from scipy.optimize import leastsq

# Optimization call with FULL Lagrangian basis

#full = 1 #lag, mesg = leastsq(residual_err, guess_full, args=(1), xtol=1.49012e-5, ftol=1.49012e-4)

# Optimization call with TRUNCATED Lagrangian basis

full = 0 lag, mesg = leastsq(residual_err2, guess_trunc, args=(0), xtol=1.49012e-5,ftol=1.49012e-5)

print mesg

''' Ux_kop = np.reshape(lag[:Nnnx_k],(Nn,nx_k),'F') Vx_kop = np.reshape(lag[Nnnx_k:(Nnnx_k + (n_t+1)nx_k)],(nx_k,n_t+1),'F') U_kop = np.reshape(lag[(Nnnx_k + (n_t+1)nx_k):(Nnnx_k + (n_t+1)nx_k + Nnn_k)],(Nn,n_k),'F') #U_kopT = U_kop.conj().T V_kop = np.reshape(lag[(Nnnx_k + (n_t+1)nx_k + Nnn_k):],(n_k,n_t+1),'F')

# Optimized truncated Lagrangian coordinate basis

X_lag = Ux_kop.dot(Vx_kop)

W_temp = U_kop.dot(V_kop)

plt.plot(x,W[:,0],'k-',x,W_temp[:,0],'g-o',x,W[:,n_t/3],'k-',x,W_temp[:,n_t/3],'g-o',x,W[:,2n_t/3],'k-',x,W_temp[:,2n_t/3],'g-o',x,W[:,n_t],'k-',x,W_temp[:,n_t],'g-o') plt.title('Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')

In [25]:
def residual_err(tilde,full):
W_interp = W.copy()
if full == 1: # FULL Lagrangian basis
Uk  = np.reshape(tilde[:Nn*n_k],(Nn,n_k),'F')
Vk  = np.reshape(tilde[Nn*n_k:],(n_k,n_t+1),'F')
## For full Lagrangian basis
X_tilde = Ux.dot(Sx_full.dot(Vx))   #Xpod_k#
elif full == 0: # TRUNCATED Lagrangian basis
Uxk = np.reshape(tilde[:Nn*nx_k],(Nn,nx_k),'F')
Vxk = np.reshape(tilde[Nn*nx_k:(Nn*nx_k + (n_t+1)*nx_k)],(nx_k,n_t+1),'F')
Uk  = np.reshape(tilde[(Nn*nx_k + (n_t+1)*nx_k):(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k)],(Nn,n_k),'F')
Vk  = np.reshape(tilde[(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k):],(n_k,n_t+1),'F')
## For truncated Lagrangian basis
X_tilde = Uxk.dot(Vxk)

W_tilde = Uk.dot(Vk)
for nn,t in enumerate(tnList[1:]):
f = interpolate.interp1d(X_tilde[:,nn], W_tilde[:,nn],fill_value="extrapolate")
W_interp[:,nn] = f(x0)
#w_interp[:,nn] = np.interp(x0, X_tilde[:,nn], w_tilde[:,nn])

return (W - W_interp).flatten('F')

In [26]:
from scipy.optimize import leastsq
# Optimization call with FULL Lagrangian basis
#full = 1
#lag, mesg = leastsq(residual_err, guess_full, args=(1), xtol=1.49012e-5, ftol=1.49012e-4)

# Optimization call with TRUNCATED Lagrangian basis
full = 0
lag, mesg = leastsq(residual_err, guess_trunc, args=(0), xtol=1.49012e-5,ftol=1.49012e-5)

print mesg

2
In [27]:
# Unwrapping output for rank-k full lagrangian basis approximation
if full == 1:
U_kop  = np.reshape(lag[:(Nn*n_k)],(Nn,n_k),'F')
#U_kopT = U_kop.conj().T
V_kop  = np.reshape(lag[(Nn*n_k):],(n_k,n_t+1),'F')
# Full Lagrangian snapshot matrix
X_lag = Ux.dot(Sx_full.dot(Vx))#X.copy()
# Unwrapping output for rank-k optimally truncated lagrangian basis approximation
elif full == 0:
Ux_kop = np.reshape(lag[:Nn*nx_k],(Nn,nx_k),'F')
Vx_kop = np.reshape(lag[Nn*nx_k:(Nn*nx_k + (n_t+1)*nx_k)],(nx_k,n_t+1),'F')
U_kop  = np.reshape(lag[(Nn*nx_k + (n_t+1)*nx_k):(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k)],(Nn,n_k),'F')
#U_kopT = U_kop.conj().T
V_kop  = np.reshape(lag[(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k):],(n_k,n_t+1),'F')
# Optimized truncated Lagrangian coordinate basis
X_lag = Ux_kop.dot(Vx_kop)
#X_full = Ux.dot(Sx_full.dot(Vx))

#print X_lag.shape

In [28]:
# Small routine to check the residual error function for bugs using the optimized bases
#guess_op = np.concatenate((Ux_kop.flatten('F'), Vx_kop.flatten('F'), U_kop.flatten('F'), V_kop.flatten('F')),axis=0)
#err = residual_err(guess_op,0)
#err = np.reshape(err,(W.shape),'F')
#plt.plot(x,err[:,n_t])

In [29]:
#Rank k approximation
Wlag = W.copy()

# optimized truncated state basis on the Lagrangian grid
W_temp = U_kop.dot(V_kop)
#w_temp = U_kop.dot(U_kop.conj().T.dot(w))
for nn,t in enumerate(tnList[1:-1]):
f = interpolate.interp1d(X_lag[:,nn+1], W_temp[:,nn+1],fill_value="extrapolate")
Wlag[:,nn+1] = f(x0)
#wlag[:,nn] = np.interp(x0, X_lag[:,nn], w_temp[:,nn])


In [30]:
plt.plot(x,W[:,0],'k-',x,Wlag[:,0],'g-o',x,W[:,n_t/3],'k-',x,Wlag[:,n_t/3],'g-o',x,W[:,2*n_t/3],'k-',x,Wlag[:,2*n_t/3],'g-o',x,W[:,n_t],'k-',x,Wlag[:,n_t],'g-o')
plt.title('Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')

Text(0.5,1,u'Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')
In [31]:
plt.plot(x,W[:,0],'k-',x,Wpod[:,0],'r-*',x,Wlag[:,0],'g-o',x,W[:,n_t/2],'k-',x,Wpod[:,n_t/2],'r-*',x,Wlag[:,n_t/2],'g-o',x,W[:,5*n_t/6],'k-',x,Wpod[:,5*n_t/6],'r-*',x,Wlag[:,5*n_t/6],'g-o')
#,x,wpod[:,n_t],'k:',x,wlag[:,n_t],'m:'
plt.title('Comparison of truth solution (black), Eulerian MOR (red stars) \n and Lagrangian MOR (green dots)')

Text(0.5,1,u'Comparison of truth solution (black), Eulerian MOR (red stars) \n and Lagrangian MOR (green dots)')
In [32]:
def plot_lag_error(index):
diff = W-Wlag
max_err=np.max(diff.flat)
min_err=np.min(diff.flat)
fig = plt.figure()
ax = fig.gca()
ax.plot(x,W[:,index]-Wlag[:,index],'k')
plt.xlabel('x')
plt.ylabel('w')
#ax.set_ylim([min_err,max_err])
plt.title('LAG error  at time $t^i=%6.3f$' % (tnList[index]))
leg=ax.legend(['$\epsilon_p^i$'],loc=1)

In [ ]:
plot_lag_error(n_t/4)


## Online evaluation - No idea what's going on here!!!

In [33]:
# Left hand side of 2nd order centrally discretized system
main_diag_online = np.ones((Nn),)
diagonals_online = [(1.0 + f2*2.0*dt/(dx**2))*main_diag_online, - f2*dt/(dx**2)*main_diag_online, - f2*dt/(dx**2)*main_diag_online]
offsets = [0, 1, -1]
#diagonals = [(1+ 2*dt*f2/dx**2)*main_diag_online, (0.5*f1*dt/dx - dt*f2/dx**2)*main_diag_online, -(0.5*f1*dt/dx + f2*dt/dx**2)*main_diag_online]
#offsets = [0, 1, -1]

# resetting coefficient matrix for dirichlet boundary
diagonals_online[0][0] = 1.0; diagonals_online[0][-1]= 1.0
diagonals_online[1][0] = 0.0
diagonals_online[2][-2]= 0.0

A_online = scipy.sparse.diags(diagonals_online, offsets, shape = (Nn,Nn), format = 'csr')

# right hand side for algebraic system
b_online = np.zeros((Nn,),'d')

In [39]:
# number of time steps
nt_online = 100
# constant time step size
dt_online =(T-t0)/float(nt_online)
DT_online = dt_online * np.ones(nt_online)

tnList_online = t0*np.ones(nt_online+1)
for i in np.arange(1,nt_online+1):
tnList_online[i] = tnList_online[i-1]+DT_online[i-1]

#U_Eul_kop = U_kop.copy()

#interpolate Lagrangian state basis from Lagrangian to Eulerian coordinates
#for nn in np.arange(n_k):
#    f = interpolate.interp1d(Ux_kop[:,nn], U_kop[:,nn],fill_value="extrapolate")
#    U_Eul_kop[:,nn] = f(x0)

Anew = A.copy()
Ahat = U_kop.conj().T.dot(A_online.dot(U_kop))
print Ahat.shape

t0 = tnList_online[0]; tn = t0; #ii=0;
Z  = np.zeros((n_k,len(tnList_online)),'d')
Wnew = np.zeros((Nn,len(tnList_online)),'d')

print U_kop.conj().T.shape
print Z.shape
print Wnew.shape

Z[:,0] = U_kop.conj().T.dot(w0)
Wnew[:,0] = w0
for nn,t in enumerate(tnList_online[1:]):
#print "step {0} solving to t={1}".format(nn+1,t)
dt=t-tn #allows for variable time step
bhat = (Z[:,nn]) #contribution from prev time step
bhat[0] = y_l # left boundary condition
bhat[-1] = y_r # right boundary condition
#print bhat.shape
bnew = (Wnew[:,nn]) #contribution from prev time step
#bnew[0] = y_l # left boundary condition
#bnew[-1] = y_r # right boundary condition
Z[:,nn+1] = scipy.sparse.linalg.spsolve(Ahat,bhat)
Wnew[:,nn+1] = scipy.sparse.linalg.spsolve(Anew,bnew)
tn=t

Znew = U_kop.dot(Z)
print Znew.shape
print X_lag.shape

Znew_E = Wnew.copy()
#interpolate Lagrangian state basis from Lagrangian to Eulerian coordinates
for nn in np.arange(nt_online):
f = interpolate.interp1d(X_lag[:,nn], Znew[:,nn],fill_value="extrapolate")
Znew_E[:,nn] = f(x0)

(3, 3) (3, 70) (3, 101) (70, 101) (70, 101) (70, 71)
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-39-eed75b7ed6af> in <module>() 53 #interpolate Lagrangian state basis from Lagrangian to Eulerian coordinates 54 for nn in np.arange(nt_online): ---> 55 f = interpolate.interp1d(X_lag[:,nn], Znew[:,nn],fill_value="extrapolate") 56 Znew_E[:,nn] = f(x0) IndexError: index 71 is out of bounds for axis 1 with size 71 
In [ ]:
plt.plot(x,Wnew[:,nt_online/5],'k-',x,Znew_E[:,nt_online/5],'g-o') #
#,x,Wnew[:,n_t/3],'k-',x,Znew[:,n_t/3],'g-o',x,Wnew[:,2*n_t/3],'k-',x,Znew[:,2*n_t/3],'g-o',x,Wnew[:,n_t],'k-',x,Znew[:,n_t],'g-o'
plt.title('Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')


## A quick demonstration of the difference in "numpy.interp" and "scipy.interpolate.interp1d" routines

In [ ]:
xx = np.arange(0,10); yy = np.array([0, 5, 1, 3, 6, 9, 2, 4, 6, 8])#10*xx+1
# grid point to extrapolate on
xxhat = xx[-1] + 10
# Using numpy interpolation with constant extension extrapolation
yyhat1 = np.interp(xxhat,xx,yy)
# Using scipy interpolation with explicit linear extrapolation
F = interpolate.interp1d(xx,yy,fill_value="extrapolate")
yyhat2 = F(xxhat)
print F(5.6)
print 'xhat = {0},yhat1 = {1}'.format(xxhat,yyhat1)
print 'xhat = {0},yhat1 = {1}'.format(xxhat,yyhat2)

# Plotting original data
plt.subplot(2, 1, 1)
plt.plot(xx,yy)
# Plotting the two types of extrapolated curves
plt.subplot(2, 1, 2)
plt.plot(np.append(xx,xxhat),np.append(yy,yyhat1),'r',np.append(xx,xxhat),np.append(yy,yyhat2),'g')