1D FEM for nonlinear diffusion coefficient

Question

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

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):
        print('iteration',str(n))
        u = uprev.copy()
        uold = u.copy() + 1
        it = 0
        while (np.linalg.norm(u-uold)>tol):
            uold=u.copy()   
            u = np.linalg.solve(odx*massmatrix(M) + np.diag(a(u))@stiffassembly(M), odx*massmatrix(M)@uprev)
            errnrm = np.linalg.norm(u-uold)
            print(errnrm)
        uprev = u.copy()
        
    
    plt.plot(x,np.r_[0,u,0],'g-o',)

---

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

Answer 1

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

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

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

Answer 2

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 $$