# 1D FEM for nonlinear diffusion coefficient

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

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

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

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

• 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$$ – Vefhug Sep 19 '20 at 7:49
• 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. – knl Sep 19 '20 at 9:40
• 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. – knl Sep 19 '20 at 9:47
• 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? – Vefhug Sep 19 '20 at 10:09
• Oh yes, sorry, you're correct. Basically, you interpolate linearly between those two values, $b_i$ and $b_{i+1}$, inside the element. – knl Sep 19 '20 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$$

• 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 – Vefhug Sep 18 '20 at 7:21
• Also, what is, explicitely, $a(u_h)$? – Vefhug Sep 18 '20 at 7:22
• You must know the function $a(u)$. So $a(u_h)$ just means evaluate it at the numerical solution. – cfdlab Sep 18 '20 at 13:09
• 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$? – Vefhug Sep 18 '20 at 13:31
• 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. – cfdlab Sep 19 '20 at 5:12