title: Matrix vector multiplication where the vector has been interpolated - Python
body: <p>I have used the finite element method to approximate the laplace equation <a href="https://i.stack.imgur.com/IA770.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/IA770.png" alt="$-u_{xx} = sin(\pi*x)$"></a> and thus have turned it into a matrix system AU = F where A is the stiffness vector and solved for U (not massively important for my question).</p>
  
  <p>I have now got my approximation U, which when i find AU i should get the vector F (or at least similar) where F is:</p>
  
  <p><a href="https://i.stack.imgur.com/QtZUA.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/QtZUA.png" alt="enter image description here"></a></p>
  
  <p>AU gives the following plot for x = 0 to x = 1 (say, for 20 nodes):</p>
  
  <p><a href="https://i.stack.imgur.com/YzoIB.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YzoIB.png" alt="enter image description here"></a></p>
  
  <p>I then need to interpolate U to a longer vector and find AU (for a bigger A too, but not interpolating that). I interpolate U by the following:</p>
  
  <pre><code>U_inter = interp1d(x,U)
  U_rich = U_inter(longer_x)
  </code></pre>
  
  <p>which seems to work okay until i multiply it with the longer A matrix:</p>
  
  <p><a href="https://i.stack.imgur.com/rm5MG.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/rm5MG.png" alt="enter image description here"></a></p>
  
  <p>It seems each spike is at a node of x (i.e. the nodes of the original U). Does anybody know what could be causing this? The following is my code to find A, U and F.</p>
  
  <pre><code>import numpy as np
  import math
  import scipy
  from scipy.sparse import diags
  import scipy.sparse.linalg
  from scipy.interpolate import interp1d
  import matplotlib
  import matplotlib.pyplot as plt
  
  def Poisson_Stiffness(x0):
      """Finds the Poisson equation stiffness matrix with any non uniform mesh x0"""
  
      x0 = np.array(x0)
      N = len(x0) - 1 # The amount of elements; x0, x1, ..., xN
  
      h = x0[1:] - x0[:-1]
  
      a = np.zeros(N+1)
      a[0] = 1 #BOUNDARY CONDITIONS
      a[1:-1] = 1/h[1:] + 1/h[:-1]
      a[-1] = 1/h[-1]
      a[N] = 1 #BOUNDARY CONDITIONS
  
      b = -1/h
      b[0] = 0 #BOUNDARY CONDITIONS
  
      c = -1/h
      c[N-1] = 0 #BOUNDARY CONDITIONS: DIRICHLET
  
      data = [a.tolist(), b.tolist(), c.tolist()]
      Positions = [0, 1, -1]
      Stiffness_Matrix = diags(data, Positions, (N+1,N+1))
  
      return Stiffness_Matrix
  
  def NodalQuadrature(x0):
      """Finds the Nodal Quadrature Approximation of sin(pi x)"""
  
      x0 = np.array(x0)
      h = x0[1:] - x0[:-1]
      N = len(x0) - 1
  
      approx = np.zeros(len(x0))
      approx[0] = 0 #BOUNDARY CONDITIONS
  
      for i in range(1,N):
          approx[i] = math.sin(math.pi*x0[i])
          approx[i] = (approx[i]*h[i-1] + approx[i]*h[i])/2
  
      approx[N] = 0 #BOUNDARY CONDITIONS
  
      return approx
  
  def Solver(x0):
  
      Stiff_Matrix = Poisson_Stiffness(x0)
  
      NodalApproximation = NodalQuadrature(x0)
      NodalApproximation[0] = 0
  
      U = scipy.sparse.linalg.spsolve(Stiff_Matrix, NodalApproximation)
  
      return U
  
  x = np.linspace(0,1,10)
  rich_x = np.linspace(0,1,50)
  U = Solver(x)
  A_rich = Poisson_Stiffness(rich_x)
  U_inter = interp1d(x,U)
  U_rich = U_inter(rich_x)
  AUrich = A_rich.dot(U_rich)
  plt.plot(rich_x,AUrich)
  plt.show()
  </code></pre>