2
$\begingroup$

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')
$\endgroup$
  • $\begingroup$ One solution of this problem (with this BC) is $u\equiv 0$. If your code is correct, this fixed -point is not attractive. Maybe there are more solutions, I don't know. The iteration is also not necessarily guaranteed to converge. Have you tried other initial values? $\endgroup$ – cos_theta Sep 22 at 16:13
  • $\begingroup$ Yes, but the situation does not change. I think that my post is correct up to the $$ -M U -B(U) U =0$$ part,right? With Mathematica, they showed that in case the equation is with $\sqrt{u}$ , the problem has a unique solution, so I edited the question changing the equation @cos_theta $\endgroup$ – Vefhug Sep 22 at 16:23
  • $\begingroup$ See also mathematica.stackexchange.com/questions/230536/… for Mathematica solution @cos_theta $\endgroup$ – Vefhug Sep 22 at 17:20
  • $\begingroup$ For $\sqrt{u}$ case, the weak form is not correct. It should be$$ \int_{-L}^L \sum_{j} u_j \phi_j(x) \phi_i'(x) \mathrm{sgn}(x)\,\mathrm{d}x + \int_{-L}^L \sqrt{ \sum_k u_k \phi_k(x) } \sum_j u_j \phi_j'(x) \phi_i'(x) \,\mathrm{d}x = 0.$$ In addition, you can try to start the iteration very close to the solution, which you know from the Mathematica.SE thread. $\endgroup$ – cos_theta Sep 22 at 17:57
  • 1
    $\begingroup$ Ah, alright. So you know $U^k$. Then what is the problem with assembly? You replace the integral by quadrature, which requires you to evaluate the integrand at individual points. You can easily do that with the square root of the sum. $\endgroup$ – Wolfgang Bangerth Sep 24 at 17:27
1
$\begingroup$

As the mathematica.se thread shows, the solution of $$ \begin{aligned}\frac{\partial}{\partial x}\left( \operatorname{sign}(x) u(x) \right) + \frac{\partial}{\partial x} \left( \sqrt{u(x)} \frac{\partial u}{\partial x}(x) \right) &= 0 & &\text{in } \Omega = (-6,6), \\ u &= 0 & &\text{on } \partial \Omega = \{-6,6\} \end{aligned}$$ is not unique. There is one non-trivial solution and the other solution is $u \equiv 0$.

Formulating the equation as $$ -\frac{\partial}{\partial x}\left( -\operatorname{sign}(x) u(x) \right) + \frac{\partial}{\partial x} \left( \sqrt{u(x)} \frac{\partial u}{\partial x}(x) \right) = 0,$$ we see that the velocity of the advection is $-\operatorname{sign}(x)$. That is, mass is always transported towards $x=0$. This also explains the shape of the solution from the mathematica.se thread, which is non-differentiable at $x=0$.

Following the usual steps, we derive the weak form $$ \lim_{a\nearrow 0} \left[ \operatorname{sign}(a)u(a)v(a) \right] - \lim_{b\searrow 0} \left[ \operatorname{sign}(b)u(b)v(b) \right] -\int_{\Omega} \operatorname{sign}(x) u(x) \frac{\partial v}{\partial x}(x)\,\mathrm{d}x - \int_{\Omega} \sqrt{u(x)} \frac{\partial u}{\partial x}(x) \frac{\partial v}{\partial x}(x) \,\mathrm{d}x= 0, $$ which simplifies to $$ -2u(0)v(0) -\int_{\Omega} \operatorname{sign}(x) u(x) \frac{\partial v}{\partial x}(x)\,\mathrm{d}x - \int_{\Omega} \sqrt{u(x)} \frac{\partial u}{\partial x}(x) \frac{\partial v}{\partial x}(x) \,\mathrm{d}x= 0 $$ provided that $u,v$ are continuous in $x=0$. Taking $u,v \in H^1_0(\Omega)$, this is indeed the case due to Sobolev embedding.

We discretize the space $H^1_0(\Omega)$ by standard hat functions $\varphi_i$ that are placed on an equidistant grid of size $h$. That is, we have $V_h = \operatorname{span}\left\{ \varphi_i : i \in \mathcal{I} \right\} \subset H^1_0(\Omega)$, where $\mathcal{I}$ is some index set.

Using this basis, we construct the matrices $A$ and $B(w)$, where $$\begin{aligned} A_{i,j} &= -2\varphi_j(0)\varphi_i(0) -\int_{\Omega} \operatorname{sign}(x) \varphi_j(x) \frac{\partial \varphi_i}{\partial x}(x)\,\mathrm{d}x, \\ B_{i,j}(w) &= - \int_{\Omega} \sqrt{w(x)} \frac{\partial \varphi_j}{\partial x}(x) \frac{\partial \varphi_i}{\partial x}(x) \,\mathrm{d}x.\end{aligned} $$ Here, the matrix $B$ still depends on some function $w \in V_h$. This gives rise to the (discrete) fixed-point problem $$A \vec{u} + B(u_h) \vec{u} = \vec{0},$$ where $\vec{u}$ denotes the coordinates of $u_h \in V_h$.

We apply a fixed-point iteration by linearizing the problem as follows:

  1. Choose $u_0 \in V_h$ and set $n = 0$.
  2. Solve $\displaystyle \left(A + B(u_n)\right) \vec{u}_{n+1} = \vec{0}$ to obtain $\vec{u}_{n+1}$.
  3. Check convergence / stopping criterion.
  4. If criterion is not satisfied, increase $n$ and go to step 2.

I've quickly hacked this scheme together in the following Python script (it is highly inefficient and doesn't even use sparse matrices). It always converges to $u \equiv 0$, even if started very close to the other solution. One can obtain a non-trivial solution if a non-zero right-hand side is applied (commented out).

#!/usr/bin/env python3

import numpy as np

def simpson(f, a,b):
    eps = np.finfo(float).eps
    # Avoid evaluating directly on the edges of the interval because of discontinuities
    return (b-a-10*eps)/6 * np.dot(np.array([1,4,1]), f(np.array([a+5*eps, (a+b)/2, b-5*eps])))

def hatFun(x, i, grid):
    if i == 0:
        center = grid[i]
        right = grid[i+1]
        return (-(x - center) / (right - center) + 1) * (x > center) * (x <= right)
    elif i == len(grid)-1:
        center = grid[i]
        left = grid[i-1]
        return (x - left) / (center-left) * (x <= center) * (x >= left)
    else:
        center = grid[i]
        left = grid[i-1]
        right = grid[i+1]
        return (x - left) / (center-left) * (x <= center) * (x >= left) + (-(x - center) / (right - center) + 1) * (x > center) * (x <= right)

def hatFunGrad(x, i, grid):
    if i == 0:
        center = grid[i]
        right = grid[i+1]
        return -1 / (right - center) * (x > center) * (x <= right)
    elif i == len(grid)-1:
        center = grid[i]
        left = grid[i-1]
        return 1 / (center-left) * (x <= center) * (x >= left)
    else:
        center = grid[i]
        left = grid[i-1]
        right = grid[i+1]
        return 1 / (center-left) * (x <= center) * (x >= left) - 1 / (right - center) * (x > center) * (x <= right)

def assembleMats(u, grid, intByParts=True):
    A = np.zeros((len(grid)-2, len(grid)-2))
    B = np.zeros((len(grid)-2, len(grid)-2))
    for i in range(1, len(grid)-1): # Test function
        idxRow = i-1
        for j in range(i-1,i+2): # Ansatz function
            if (j == 0) or (j == len(grid)-1):
                # Early out for non-overlapping support
                continue
            idxCol = j-1

            if intByParts:
                if ((grid[i-1] < 0) and (grid[i+1] <= 0)):
                    A[idxRow, idxCol] += simpson(lambda x: hatFun(x, j, grid) * hatFunGrad(x, i, grid), grid[i-1], grid[i])
                    A[idxRow, idxCol] += simpson(lambda x: hatFun(x, j, grid) * hatFunGrad(x, i, grid), grid[i], grid[i+1])
                elif (grid[i-1] >= 0):
                    A[idxRow, idxCol] -= simpson(lambda x: hatFun(x, j, grid) * hatFunGrad(x, i, grid), grid[i-1], grid[i])
                    A[idxRow, idxCol] -= simpson(lambda x: hatFun(x, j, grid) * hatFunGrad(x, i, grid), grid[i], grid[i+1])
                else: # grid[i-1] < 0, grid[i] == 0, grid[i+1] > 0

                    # \int_{-h}^{0} d/dx( sign(x) phi_j ) * phi_i dx
                    #   = [sign * phi_j * phi_i]_{-h}^{0} - \int_{-h}^{0} sign(x) phi_j * dphi_i/dx dx
                    #   = [-phi_j * phi_i]_{-h}^{0} + \int_{-h}^{0} phi_j * dphi_i/dx dx
                    #   = -phi_j(0) * phi_i(0) + \int_{-h}^{0} phi_j * dphi_i/dx dx
                    A[idxRow, idxCol] += simpson(lambda x: hatFun(x, j, grid) * hatFunGrad(x, i, grid), grid[i-1], grid[i]) \
                        -hatFun(0, j, grid) * hatFun(0, i, grid)

                    # \int_{0}^{h} d/dx( sign(x) phi_j ) * phi_i dx
                    #   = [sign * phi_j * phi_i]_{0}^{h} - \int_{0}^{h} sign(x) phi_j * dphi_i/dx dx
                    #   = [phi_j * phi_i]_{0}^{h} - \int_{0}^{h} phi_j * dphi_i/dx dx
                    #   = -phi_j(0) * phi_i(0) - \int_{0}^{h} phi_j * dphi_i/dx dx
                    A[idxRow, idxCol] += -simpson(lambda x: hatFun(x, j, grid) * hatFunGrad(x, i, grid), grid[i], grid[i+1]) \
                        -hatFun(0, j, grid) * hatFun(0, i, grid)
            else:
                if ((grid[i-1] < 0) and (grid[i+1] <= 0)):
                    A[idxRow, idxCol] -= simpson(lambda x: hatFunGrad(x, j, grid) * hatFun(x, i, grid), grid[i-1], grid[i])
                    A[idxRow, idxCol] -= simpson(lambda x: hatFunGrad(x, j, grid) * hatFun(x, i, grid), grid[i], grid[i+1])
                elif (grid[i-1] >= 0):
                    A[idxRow, idxCol] += simpson(lambda x: hatFunGrad(x, j, grid) * hatFun(x, i, grid), grid[i-1], grid[i])
                    A[idxRow, idxCol] += simpson(lambda x: hatFunGrad(x, j, grid) * hatFun(x, i, grid), grid[i], grid[i+1])
                else: # grid[i-1] < 0, grid[i] == 0, grid[i+1] > 0
                    A[idxRow, idxCol] -= simpson(lambda x: hatFunGrad(x, j, grid) * hatFun(x, i, grid), grid[i-1], grid[i])
                    A[idxRow, idxCol] += simpson(lambda x: hatFunGrad(x, j, grid) * hatFun(x, i, grid), grid[i], grid[i+1])

            B[idxRow, idxCol] = simpson(lambda x: np.sqrt( u[i-1] * hatFun(x, i-1, grid) + u[i] * hatFun(x, i, grid) + u[i+1] * hatFun(x, i+1, grid) ) * hatFunGrad(x, i, grid) * hatFunGrad(x, j, grid), grid[i-1], grid[i]) \
                + simpson(lambda x: np.sqrt( u[i-1] * hatFun(x, i-1, grid) + u[i] * hatFun(x, i, grid) + u[i+1] * hatFun(x, i+1, grid) ) * hatFunGrad(x, i, grid) * hatFunGrad(x, j, grid), grid[i], grid[i+1])

    return (A, -B)

def assembleVec(grid, f):
    v = np.zeros((len(grid)-2,))
    for i in range(1, len(grid)-1):
        idxRow = i-1
        v[idxRow] = simpson(lambda x: f(x) * hatFun(x, i, grid), grid[i-1], grid[i])
        v[idxRow] += simpson(lambda x: f(x) * hatFun(x, i, grid), grid[i], grid[i+1])

    return v

def fixedPoint(u0, rhs, grid, intByParts=False):
    nFixPoint = 50
    tol = 1e-10
    for i in range(nFixPoint):
        (A,B) = assembleMats(u0, grid, intByParts=intByParts)

        res = np.dot(A, u0[1:-1]) + np.dot(B, u0[1:-1]) - rhs
        resSq = np.sqrt(np.dot(res,res))
        print('Iter {:2d}: Residual: {:e}'.format(i, resSq))

        if resSq <= tol:
            break

        # Solve inner nodes
        un = np.linalg.solve(A+B, rhs)
        # Add outer nodes (Dirichlet BCs)
        u0 = np.r_[0, un, 0]
    return u0


# Number of points has to be odd (we need 0.0 as grid point)
grid = np.linspace(-6, 6, 11)

# Interpolation of true solution at nodal points
#u0 = np.array([0.0, 0.3600, 1.440, 3.240, 5.760, 9.000, 5.760, 3.240, 1.440, 0.3600, 0.0])

# L2 projection of solution to finite dimensional space
#u0 = np.array([0.0, 0.5040, 1.800, 3.960, 6.984, 9.432, 6.984, 3.960, 1.800, 0.5040, 0.0])

u0 = np.ones(len(grid),)

# Enforce Dirichlet BCs for initial guess
u0[0] = 0.0
u0[-1] = 0.0

# Select right hand side
rhs = np.zeros((len(grid)-2,))
#rhs = assembleVec(grid, lambda x: -np.sqrt(x + 6))

u = fixedPoint(u0, rhs, grid, intByParts=False)
uIBP = fixedPoint(u0, rhs, grid, intByParts=True)

import matplotlib.pyplot as plt
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.set_title('Solution')
ax1.plot(grid,u)
ax1.plot(grid,uIBP)
ax1.legend(['W/o IntByParts', 'W/ IntByParts'])

ax2 = fig.add_subplot(212)
ax2.set_title('Difference of solutions')
ax2.plot(grid,u-uIBP)

plt.show()

plt.plot(grid,u)
plt.show()

I'd suggest pseudo time stepping (or pseudo-transient continuation) started from a non-zero initial guess in order to compute the other non-trivial solution.

Here's why (please correct me if I'm wrong): Considering the solution as the steady state of the time-dependent equation, we see that the diffusive term (distribution of mass) exactly balances the advective term (transport towards $x=0$). Thus, in the steady state, no mass can enter or exit the system due to the boundary conditions and the flow field. In the transient phase, mass can still enter or exit the system as needed in order to reach the steady state. Therefore, a time stepping based method seems more appropriate to me than fixed-point or some kind of Newton's iteration.

For the fixed-point iteration, I suspect that $A + B(w)$ is always invertible, except for $w \in H^1_0$ being the non-trivial solution. Since we cannot exactly represent this non-trivial solution in $V_h$, we always end up with $u \equiv 0$. Thus, fixed-point iteration is not suitable here.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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