I need to solve with linear finite elements the equation $$\frac{\partial }{\partial x}\Bigl(\text{sgn}(x) u \Big) +\frac{\partial}{\partial x} \Bigl[ \sqrt{u} \frac{\partial u}{\partial x} \Bigr] =0$$
with boundary conditions $u(-L)=u(L)=0$ where $L=6$
(It's the steady state version of the equation described here: Non-Linear advection diffusion with nondifferetiable advection term)
I take $v \in H_0^1(-L,L)$ and after the usual steps I obtain $$- \int u(x) \text{sgn}(x) v' dx - \int \sqrt{u} u' v' dx = 0$$
Then, using linear finite elements: $$- \int \sum_{j} u_j \phi_j(x) \phi_i'(x) \text{sgn}(x)dx - \int \Bigl( \sum_k \sqrt{u_k \phi_k} \Bigr) \sum_j u_j \phi_j' \phi_i' dx = 0$$
which leads to the non-linear system (setting $U=[u_0,\ldots,u_N]$)
$$-C U -A(U) U =$$
where $(C)_{ij} = \int \sum_{j} u_j \phi_j(x) \phi_i'(x) \text{sgn}(x)dx $
and $\Bigl(A(U)\Bigr)_{ij} =\int \Bigl( \sum_k \sqrt{u_k \phi_k} \Bigr) \sum_j u_j \phi_j' \phi_i' dx$
Now, I want to solve this non-linear equation with fix-point iterations, so I set $$CU^{k+1} = -A(U^k)U^k$$ and solve iteratively those linear systems.
The problem: unfortunately, the fix-point iteration gives me NaN and I can't find the solution. Is it because the problem is ill-posed, or did I do something wrong with my idea of fixpoint iterations?
After @cos_theta's comment, I modified my code with the right weak formulation, but still the solution can't be found. Basically, I made two functions, one where I assemble the matrix $A(U)$, and the other one where I assemble the matrix $C$. Then I have the fixed-point iteration loop.
In particular, the matrix $A(U)$ corresponds to $$\int \sqrt{ \sum_k u_k \phi_k } \sum_j u_j \phi_j' \phi_i' dx = 0$$
so it's tridiagonal and, for instance, the diagonal entry is $$\int_{x_{i-1}}^{x_i} \sqrt{u_{i-1}}\sqrt{\phi_{i-1}} \frac{1}{h^2}dx + \int_{x_i}^{x_{i+1}} \sqrt{u_{i+1}} \sqrt{\phi_{i+1}} \frac{1}{h^2}dx + \int_{x_{i-1}}^{x_{i+1}} \sqrt{u_i} \sqrt{\phi_i} \frac{1}{h^2} dx $$
where the values $\sqrt{u_{i-1}}$, $\sqrt{u_i}$, $\sqrt{u_{i+1}}$ are given by the previous iteration.
For the matrix $C$, I have that $$C_{ii}= \int_{x_{i-1}}^{x_i} \frac{1}{h} \phi_i \text{sgn}(x) dx + \int_{x_i}^{x_{i+1}} \frac{-1}{h} \phi_i \text{sgn}(x)dx$$ If the interval does not contain $x=0$, then $C_{ii}=0$. Otherwise, as shown in the linked answer, the entry that contains $x=0$ is $-1$. So the resulting matrix is like this
$$C = \begin{pmatrix}0 & \frac{1}{2} & 0 & 0 & 0 \\ -\frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 \\ 0 & -\frac{1}{2} & -1 & -\frac{1}{2} & 0 \\ 0 & 0 & \frac{1}{2} & 0 & -\frac{1}{2} \\ 0 & 0 & 0 & \frac{1}{2} & 0\end{pmatrix}$$
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
L = 6
def stiffassembly(a,M):
# a is the vector containg the previous solution. It's long M+1, it takes also boundary values in order to assemble the matrix
x = np.linspace(-L,L,M+1)
diag = np.zeros(M-1) #x_1,...,x_M-1 (M-1)
supr = np.zeros(M-2)
h = x[1]-x[0]
c = 1.0/(h**2)
for i in range(1,M):
diag[i-1] = a[i-1]*c*integrate.quad(lambda t: np.sqrt((x[i]-t)/h),x[i-1],x[i])[0] + a[i+1]*c*integrate.quad(lambda t: np.sqrt((t-x[i])/h),x[i],x[i+1])[0] + a[i]*( integrate.quad(lambda t: np.sqrt((t-x[i-1])/h),x[i-1],x[i])[0] + integrate.quad(lambda t: np.sqrt((x[i+1]-t)/h),x[i],x[i+1])[0] )
for k in range(1,M-1):
supr[k-1] = a[k]*(-c)*integrate.quad(lambda t:np.sqrt((x[k+1]-t)/h),x[k],x[k+1])[0] + a[k+1]*(-c)*integrate.quad(lambda t: np.sqrt((t - x[k])/h),x[k],x[k+1])[0]
A = np.diag(supr,-1) + np.diag(diag,0) + np.diag(supr,+1)
return A
def Cmatrix(M):
x = np.linspace(-L,L,M+1)
diag = np.zeros(M-1)
subd = np.zeros(M-2)
supr = np.zeros(M-2)
h = x[1]-x[0]
c = 1.0/(h**2)
for i in range(1,M):
diag[i-1] = c*integrate.quad(lambda t: np.sign(t)*(t-x[i-1]),x[i-1],x[i])[0] - c*integrate.quad(lambda t: np.sign(t)*(x[i+1] - t),x[i],x[i+1])[0]
for k in range(1,M-1):
supr[k-1] = c*integrate.quad(lambda t:np.sign(t)*(x[k+1]-t),x[k],x[k+1])[0]
subd[k-1] = -c*integrate.quad(lambda t: np.sign(t)*(t-x[k]),x[k],x[k+1])[0]
C = np.diag(supr,-1) + np.diag(diag,0) + np.diag(subd,+1)
return C
a = lambda w: np.real(np.sqrt(w))
M = 100
x = np.linspace(-L,L,M+1)
tol = 1e-14
ts = 1000
bc = np.array([0,0])
uold = np.ones(M-1)
it = 0
errnrm = 1
C = Cmatrix(M)
while (errnrm>tol):
it+=1
u = np.linalg.solve(C,-stiffassembly(a(np.r_[bc[0],uold,bc[1]]), M)@uold)
errnrm = np.linalg.norm(u-uold)
uold = u.copy()
print(errnrm)
plt.figure()
plt.plot(x,np.r_[bc[0],u,bc[1]],'-')
plt.xlabel('x')
