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. 935
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.
```python 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. 935