Beginner Finite Elemente Code does not solve equation properly

Question

I am trying to write the code for solving the extremely difficult differential equation: x' = 1 with the finite element method. As far as I understood, I can obtain the solution u as

with the basis functions phi_i(x), while I can obtain the u_i as the solution of the system of linear equations:

with the differential operator D (here only the first derivative). As a basis I am using the tent function:

    def tent(l, r, x):
    m = (l + r) / 2 
    if x >= l and x <= m:
        return (x - l) / (m - l)
    elif x < r and x > m:
        return (r - x) / (r - m)
    else:
        return 0

def tent_half_down(l,r,x):
    if x >= l and x <= r:
        return (r - x) / (r - l)
    else:
        return 0

def tent_half_up(l,r,x):
    if x >= l and x <= r:
        return (x - l) / (r - l)
    else:
        return 0

def tent_prime(l, r, x):
    m = (l + r) / 2 
    if x >= l and x <= m:
        return 1 / (m - l)
    elif x < r and x > m:
        return 1 / (m - r)
    else:
        return 0

def tent_half_prime_down(l,r,x):
    if x >= l and x <= r:
        return - 1 / (r - l)
    else:
        return 0

def tent_half_prime_up(l, r, x):
    if x >= l and x <= r:
        return 1 / (r - l)
    else:
        return 0

def sources(x):
    return 1

Discretizing my space:

n_vertex = 30
n_points = (n_vertex-1) * 40
space = (0,5)
x_space = np.linspace(space[0],space[1],n_points)
vertx_list = np.linspace(space[0],space[1], n_vertex)
tent_list = np.zeros((n_vertex, n_points))
tent_prime_list = np.zeros((n_vertex, n_points))

tent_list[0,:] = [tent_half_down(vertx_list[0],vertx_list[1],x) for x in x_space]
tent_list[-1,:] = [tent_half_up(vertx_list[-2],vertx_list[-1],x) for x in x_space]
tent_prime_list[0,:] = [tent_half_prime_down(vertx_list[0],vertx_list[1],x) for x in x_space]
tent_prime_list[-1,:] = [tent_half_prime_up(vertx_list[-2],vertx_list[-1],x) for x in x_space]
for i in range(1,n_vertex-1):
    tent_list[i, :] = [tent(vertx_list[i-1],vertx_list[i+1],x) for x in x_space]
    tent_prime_list[i, :] = [tent_prime(vertx_list[i-1],vertx_list[i+1],x) for x in x_space]

Calculating the system of linear equations:

b = np.zeros((n_vertex))
A = np.zeros((n_vertex,n_vertex))
for i in range(n_vertex):
    b[i] = np.trapz(tent_list[i,:]*sources(x_space))
    for j in range(n_vertex):
        A[j, i] = np.trapz(tent_prime_list[j] * tent_list[i])

And then solving and reconstructing it

u = np.linalg.solve(A,b)
sol = tent_list.T.dot(u)

But it does not work, I am only getting some up and down pattern. What am I doing wrong?

Answer 1

First, a couple of comments on terminology and notation:

1) You are using the weak formulation, though you've done this implicitly. A formulation being "weak" has nothing to do with the order of derivatives involved. It is weak because you are not satisfying the differential equation exactly at every location. FE minimizes the weighted residual of the solution, integrated over the domain. The functions phi_j actually discretize the weighting function. The difference when you only have first-order derivatives is that you don't have to apply the Gauss divergence theorem (which simplifies to integration by parts for one dimension) to eliminate second-order derivatives. You can tell this wasn't done because phi_j is not differentiated in the LHS.

2) I would suggest not using "A" as the differential operator. You also use this symbol for the global system matrix, so your notation is inconsistent. People often use "D", since this fits better to the idea that it is used for differentiation.

Secondly, about your implementation:

3) You are using way more integration points than necessary. Your elements use linear interpolation functions, which means you only need one integration point located at the center of the element to evaluate the integral exactly. Look into the details of Gauss quadrature to see why. Also, you've specified the number of integration points as a multiple of the number of nodes. This should be done as a multiple of the number of elements instead (in your case, n_vertex-1), because the elements are the domains on which you're integrating.

4) You have built your system by simply removing the two end nodes from the formulation. This isn't the correct way to specify boundary conditions. I would suggesting building the full system first and using one of the typical methods for applying Dirichlet boundary conditions. Also, think about what constraining two nodes would imply for the differential equation you're trying to solve. What function exists that satisfies x' = 1, x(0) = 0, x(5) = 0? You have overconstrained the system by trying to apply 2 boundary conditions to a first-order differential equation.

Unfortunately, there isn't a small tweak that can be made to get the code to work, but I hope the comments above help you rethink your approach.

EDIT to address your changes:

1) Assuming the matrix A is addressed with A[row,col], then your indices are backwards. You should be integrating with A[i,j] = ...

2) A simple way to apply a constraint is to replace one row with the constraint desired. If you want x(0) = 0, for example, set A[0,j] = 0 for all j, then set A[0,0] = 1 and set b[0] = 0. This substitutes one of the equations with u_0 = 0. Do this after integrating.