I want to solve with linear finite elements the equation $$\partial_t u = \partial_{x}(a(u)\partial_xu)$$ in the domain $t \in [0,1]$ and $x \in [-L,L]$. Here $a(u)$ is just a function of $u$.

Applying the weak formulation with $u(t,x)=\sum_{j} u_j(t) \varphi_j(x)$, I obtain $$\partial_t u_j(t) \int_{-L}^{L}\varphi_i(x)\varphi_j(x)dx = - \int_{-L}^{L} a\Bigl( \sum_j u_j(t) \varphi_j(x) \Bigr) \Bigl( \sum_k u_k(t) \varphi_k^{'}(x) \Bigr) \varphi_{i}^{'}(x)dx$$

  • The l.h.s is no problem because it is $M \dot{U}(t)$, where $(M)_{ij}=\int_{-L}^{L} \varphi_i(x) \varphi_j(x)dx$ and $U(t)=[u(x_1,t),\ldots,u(x_N,t)]^{T}$

  • My big problem is on the r.h.s. I don't know how to handle that double summation so that I have a function of $U(t)$, because I obtain a tensor $B_{ijk}=\int_{-L}^{L} \varphi_i \varphi_j \varphi_k^{'}$ ( there has already been a question about this) but I can't understand how to solve this in practice on a computer.

As described in the linked question, I will obtain $$M \dot{U} = (BU)U$$ but this seems just formal to me, because of that tensor. Any help is highly appreciated

EDIT after knl answer:

@knl I have a question about the root-finding step:

After time discretization, I have $u_n(x)$, therefore the problem is still continuous in space. From the the scalar prodcut $$ (\delta^{-1} u_{k,n}, v) + (a(u_{k-1,n}) \partial_x u_{k,n}, \partial_x v) = (\delta^{-1}u_{n-1}, v) $$ I want to find how to compute the solution by fix point iteration

Let $A$ the usual "stiffness matrix" and $M$ the "mass matrix":

$$\delta^{-1} M u_k^n + a(u_{k-1}^n) A u_k^n = \delta^{-1}M u^{n-1}$$ where $u_k^{n}$ is the coefficients vector and $k$ is the index referring to the fix-point iteration.

Therefore, I iteratively find $u_k^n$ by solving the linear systems $$(\delta^{-1} M + a(u_{k-1}^n) A)u_k^n = \delta^{-1} M u^{n-1}$$

What I obtain after integration up to time $t=1$ is enter image description here

which is slightly different from yours. I can't understand if there's an error in my code, because the fixed point iterations seem to work.

    import numpy as np
    import matplotlib.pyplot as plt
    def stiffassembly(M):
        x = np.linspace(0,1,M+1)
        diag = np.zeros(M-1) #x_1,...,x_M-1 (M-1)
        subd = np.zeros(M-2) 
        supr = np.zeros(M-2)
        h = np.diff(x)
        for i in range(1,M):
            diag[i-1] = 1/h[i-1] +1/h[i]
        for k in range(1,M-1):
            supr[k-1] = -1/h[k]
            subd[k-1] = -1/h[k]
        A = np.diag(subd,-1) + np.diag(diag,0) + np.diag(supr,+1)
        return A
    def massmatrix(N):
        x = np.linspace(0,1,N+1)
        diag = np.zeros(N-1) #x_1,...,x_M-1 (M-1)
        subd = np.zeros(N-2) 
        supr = np.zeros(N-2)
        h = np.diff(x)
        for i in range(1,N):
            diag[i-1] = (h[i-1] + h[i])/3
        for k in range(1,N-1):
            supr[k-1] = h[k]/6
            subd[k-1] = h[k-1]/6
        M = np.diag(subd,-1) + np.diag(diag,0) + np.diag(supr,+1)
        return M
    def inidata(x):
        return np.sin(np.pi*x)
    a = lambda w: (1. * w) ** 2
    M = 50
    x = np.linspace(0,1,M+1)
    delta = 0.001
    odx = 1.0/delta
    tol = 1e-14
    uprev = inidata(x[1:-1])
    ts = 1000 #integration up to t=1.0
    for n in range(ts):
        u = uprev.copy()
        uold = u.copy() + 1
        it = 0
        while (np.linalg.norm(u-uold)>tol):
            u = np.linalg.solve(odx*massmatrix(M) + np.diag(a(u))@stiffassembly(M), odx*massmatrix(M)@uprev)
            errnrm = np.linalg.norm(u-uold)
        uprev = u.copy()

EDIT (last one) Actually what I have before is wrong, before I do not update the matrix in the fix point iteration. Fixing this, i.e. changing the function stiffassembly, i obtain the following, which is right, as it has been "three"-checked (with @knl fem solver, with a finite difference approach, and with mathematica):

enter image description here

  • $\begingroup$ @knl I don't know how to handle that tensor for a practical implementation $\endgroup$
    – Vefhug
    Commented Sep 18, 2020 at 8:12
  • $\begingroup$ I think that your edit should be an answer. $\endgroup$
    – nicoguaro
    Commented Sep 19, 2020 at 14:34

2 Answers 2


You would need to linearize the problem. I prefer to do it before discretization but it's possible to do also after discretization. (I'm a bit skeptical of linearization after discretization because I have never looked into the details. In general, discretization and linearization steps do not commute.)

In the following I assume that the equation is actually $\partial_t u = \partial_x(a(u) \partial_x u)$ and that you have the boundary condition $u=0$.

The weak form is $$(\partial_t u, v) = -(a(u) \partial_x u, \partial_x v).$$ I prefer to first do the time discretization so that you see the structure of the resulting problem. E.g., implicit Euler method leads to $$(\delta^{-1}(u_n - u_{n-1}), v) = -(a(u_n) \partial_x u_n, \partial_x v),$$ or, equivalently, $$(\delta^{-1} u_n, v) + (a(u_n) \partial_x u_n, \partial_x v) = (\delta^{-1}u_{n-1}, v),$$ where $n$ runs over the time steps and $\delta > 0$ is the size of the step. The equation is still nonlinear in $u_n$ and you must linearize. One option is to do a fixed-point iteration (inside each time step $n$) by repeatedly finding $u_{k,n}$ from $$(\delta^{-1} u_{k,n}, v) + (a(u_{k-1,n}) \partial_x u_{k,n}, \partial_x v) = (\delta^{-1}u_{n-1}, v),$$ where $k$ runs over the linearization steps and $u_{k-1,n}$ is the function from the previous iteration. Notice how you now have two iterations: one for time discretization and one for linearization.

I made an example case with $u(x) = \sin(\pi x)$ and solved it using the code I know the best (i.e. my own, you can install it in Python using pip install scikit-fem==2.0.0 if you want to run it):

from skfem import *
from skfem.helpers import *
from skfem.visuals.matplotlib import *
import numpy as np

m = MeshLine(); m.refine(5)
basis = InteriorBasis(m, ElementLineP2())
a = lambda w: (1. * w) ** 2
bilinf_stiffness = BilinearForm(lambda u, v, w: a(w['u_prev']) * dot(grad(u), grad(v)))
delta = 0.01
M = BilinearForm(lambda u, v, w: 1. / delta * u * v).assemble(basis)
load = LinearForm(lambda v, w: 1. / delta * w['u_prev'] * v)

u = project(lambda x: np.sin(np.pi * x[0]), basis_to=basis)
plot(basis, u)
for n in range(100): # 100 time steps 
    b = load.assemble(basis, u_prev=basis.interpolate(u))
    for k in range(250): # 250 linearization loops
        A = bilinf_stiffness.assemble(basis, u_prev=basis.interpolate(u))
        u = solve(*condense(A + M, b, D=m.boundary_nodes()))
    print("iteration {}".format(n))
plot(basis, u)

This gives the following two pictures (initial condition and the result at $t=1$): Initial field Solution at t=1

There are obviously lots of alternative ways of doing this, but this should give you the general idea.

  • $\begingroup$ Now it's clear. Just to conclude, I've never done such an assembly. So, basically, when I evaluate at integration points, $a(u_{k-1,n}) $ is a vector with components, say $\vec{a}$. Then, the integral becomes $$\int_0^1 \sum_j a_j \phi_j \sum_j u_j \phi _j^{'} \phi_i^{i} dx$$ right? Or is is just $$\int_0^1 \sum_j a_j \phi_j \phi _j^{'} \phi_i^{'} dx$$ $\endgroup$
    – Vefhug
    Commented Sep 19, 2020 at 7:49
  • $\begingroup$ Looking at your code, it seems that you are not using integration points at all. Instead you have already evaluated the integrals analytically with the assumption that $A_{ij} = \int_0^1 \phi_j^\prime(x) \phi_i^\prime(x) dx$. Now that you have $A_{ij} = \int_0^1 b(x) \phi_j^\prime(x) \phi_i^\prime(x) dx$ (where $b(x) = a(u_{k-1,n}(x))$) this becomes impossible unless you make further assumptions, e.g., set $b(x)$ to the mean of $a(u_{k-1,n}(x))$ over each element. This seems to be a common way to obtain analytical integrals, and usually works quite well if $b(x)$ is smooth enough. $\endgroup$
    – knl
    Commented Sep 19, 2020 at 9:40
  • $\begingroup$ Let me add that I think it might be possible to integrate analytically also when $b(x)$ is piecewise-linear. Denote by $b_i$ the value of $b(x)$ at $i$'th node $x_i$ which corresponds to $\phi_i$. Then you would have, e.g., $\int_0^1 \left(\frac{x-x_i}{x_{i+1}-x_i}b_i + \frac{x_i-x}{x_{i+1}-x_i} b_{i+1}\right) \phi_i^\prime \phi_{i+1}^\prime dx$ for the off-diagonals which I think can be integrated analytically. Nevertheless, most of the generic finite element codes are usually implemented with the help of numerical integration. $\endgroup$
    – knl
    Commented Sep 19, 2020 at 9:47
  • $\begingroup$ thanks, now I understood how to move. By the wayI was writing explicitely the integrals for the assembly, I think that your integral should be $$\int_0^1 \left(\frac{x_{i+1}-x}{x_{i+1}-x_i}b_i + \frac{x-x_i}{x_{i+1}-x_i} b_{i+1}\right) \phi_i^\prime \phi_{i+1}^\prime dx$$, right? $\endgroup$
    – Vefhug
    Commented Sep 19, 2020 at 10:09
  • $\begingroup$ Oh yes, sorry, you're correct. Basically, you interpolate linearly between those two values, $b_i$ and $b_{i+1}$, inside the element. $\endgroup$
    – knl
    Commented Sep 19, 2020 at 10:13

The equation you have is $$ M \dot U = B(U) U $$ where $B$ is a matrix which depends on the solution and is given by $$ B_{ij} = \int a(u_h) \phi_i' \phi_j' dx $$

  • $\begingroup$ could you expand a bit your answer? It's not clear to me how you got the $B_{ij}$, i.e. a matrix and not a tensor, as Wolfgang wrote in the linked question $\endgroup$
    – Vefhug
    Commented Sep 18, 2020 at 7:21
  • $\begingroup$ Also, what is, explicitely, $a(u_h)$? $\endgroup$
    – Vefhug
    Commented Sep 18, 2020 at 7:22
  • $\begingroup$ You must know the function $a(u)$. So $a(u_h)$ just means evaluate it at the numerical solution. $\endgroup$
    – cfdlab
    Commented Sep 18, 2020 at 13:09
  • $\begingroup$ My bad, sorry! I have a question about that term. When I discretize, $a(u_h)$ is a vector. So I have, after calling $A$ the matrix with $(A)_{ij} = \int \phi_i' \phi_j' dx$ $$\vec{a} A$$ where the multiplication is entrywise. Should this be written in a program as $\text{diag(a(uh)) * A$? $\endgroup$
    – Vefhug
    Commented Sep 18, 2020 at 13:31
  • $\begingroup$ Its a nonlinear problem. You have to integrate with $a(u_h)$ inside the integral. The matrix depends on the solution, so you have to recompute it every time $u_h$ changes. $\endgroup$
    – cfdlab
    Commented Sep 19, 2020 at 5:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.