Matrix vector multiplication where the vector has been interpolated - Python

Question

I have used the finite element method to approximate the laplace equation 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).

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:

AU gives the following plot for x = 0 to x = 1 (say, for 20 nodes):

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:

U_inter = interp1d(x,U)
U_rich = U_inter(longer_x)

which seems to work okay until i multiply it with the longer A matrix:

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.

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

Answer 1

comment 1:

I added a Stiffness_Matrix = Stiffness_Matrix.tocsr() statement to avoid an efficiency warning. FE calculations are complex enough that I'll have to print out some intermediate values before I can identify what is going on.

comment 2:

plt.plot(rich_x,A_rich.dot(Solver(rich_x))) plots nice. The noise you get is the result of the difference between the inperpolated U_rich and the true solution: U_rich-Solver(rich_x).

comment 3:

I don't think there's a problem with your code. The problem is with idea that you can test an interpolation this way. I'm rusty on FE theory, but I think you need to use the shape functions to interpolate, not a simple linear one.

comment 4:

Intuitively, with A_rich.dot(U_rich) you are asking, what kind of forcing F would produce U_rich. Compared to Solver(rich_x), U_rich has flat spots, regions where it's value is less than the true solution. What F would produce that? One that is spiky, with NodalQuadrature(x) at the x points, but near zero values in between. That's what your plot is showing.

A higher order interpolation will eliminate the flat spots, and produce a smoother back calculated F. But you really need to revisit the FE theory.

You might find it instructive to look at

plt.plot(x,NodalQuadrature(x))
plt.plot(rich_x, NodalQuadrature(rich_x))

The second plot is much smoother, but only about 1/5 as high.

Better yet look at:

plt.plot(rich_x,AUrich,'-*')  # the spikes
plt.plot(x,NodalQuadrature(x),'o')  # original forcing
plt.plot(rich_x, NodalQuadrature(rich_x),'+') # new forcing

In the model the forcing isn't continuous, it is a value at each node. With more nodes (rich_x) the magnitude at each node is less.