SharedScalar_Convection_Diffusion.ipynbOpen in CoCalc
Author: Corey Trahan
Views : 15

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

%%latex

where f1=1f_1 = 1, f2=1/Pef_2 = 1/Pe, w(y,0)=0.5e(y0.3)2/0.052w(y,0) = 0.5e^{-(y-0.3)^2/0.05^2}, w(0,t)=0w(0,t) = 0.

Case 1: Pe=Pe = \infty (scalar convection) and Case 2: Pe=103Pe = 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

S\mathbf{{S}} is a diagonal matrix containing the eigenvalues and U\mathbf{{U}} is a square matrix containing the eigenvectors of the kernel Kt\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

S^\mathbf{\hat{S}} is a diagonal matrix containing the eigenvalues and U^\mathbf{\hat{U}} is a square matrix containing the eigenvectors of the kernel Kxt\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

NO USE: Computing a non-optimized approximation by interpolating the reduced Eulerian state basis on the reduced Lagrangian coordinate basis

NO USE: Computing a non-optimized approximation by interpolating the reduced Eulerian state basis on the FULL Lagrangian coordinate basis

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')