Finite Elemente (1D) für ein nichtlineares stationäres Problem

Question

Ich muss die Gleichung mit linearen finiten Elementen lösen $$\frac{\partial }{\partial x}\Bigl(\text{sgn}(x) u \Big) +\frac{\partial}{\partial x} \Bigl[ \sqrt{u} \frac{\partial u}{\partial x} \Bigr] =0$$

mit Randbedingungen $u(-L)=u(L)=0$ wo $L=6$

(Es ist die stationäre Version der hier beschriebenen Gleichung: Nichtlineare Advektionsdiffusion mit nicht differenzierbarem Advektionsterm )

---

ich nehme $v \in H_0^1(-L,L)$ und nach den üblichen Schritten erhalte ich $$- \int u(x) \text{sgn}(x) v' dx - \int \sqrt{u} u' v' dx = 0$$

Dann mit linearen finiten Elementen: $$- \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$$

was zum nichtlinearen System führt (Einstellung$U=[u_0,\ldots,u_N]$)

$$-C U -A(U) U =$$

wo $(C)_{ij} = \int \sum_{j} u_j \phi_j(x) \phi_i'(x) \text{sgn}(x)dx $

und $\Bigl(A(U)\Bigr)_{ij} =\int \Bigl( \sum_k \sqrt{u_k \phi_k} \Bigr) \sum_j u_j \phi_j' \phi_i' dx$

Jetzt möchte ich diese nichtlineare Gleichung mit Fixpunktiterationen lösen , also setze ich$$CU^{k+1} = -A(U^k)U^k$$ und iterativ diese linearen Systeme lösen.

Das Problem: Leider gibt mir die Fixpunkt-Iteration NaNund ich kann keine Lösung finden. Liegt es daran, dass das Problem schlecht gestellt ist, oder habe ich mit meiner Idee von Fixpunkt-Iterationen etwas falsch gemacht?

---

Nach dem Kommentar von @ cos_theta habe ich meinen Code mit der richtigen schwachen Formulierung geändert, aber die Lösung kann immer noch nicht gefunden werden. Grundsätzlich habe ich zwei Funktionen erstellt, eine, bei der ich die Matrix zusammensetze$A(U)$und der andere, wo ich die Matrix zusammenbaue $C$. Dann habe ich die Festpunkt-Iterationsschleife.

Insbesondere die Matrix $A(U)$ entspricht $$\int \sqrt{ \sum_k u_k \phi_k } \sum_j u_j \phi_j' \phi_i' dx = 0$$

so ist es tridiagonal und zum Beispiel ist der diagonale Eintrag $$\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 $$

wo die Werte $\sqrt{u_{i-1}}$, $\sqrt{u_i}$, $\sqrt{u_{i+1}}$ sind durch die vorherige Iteration gegeben.

Für die Matrix $C$, Ich habe das $$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$$ Wenn das Intervall nicht enthält $x=0$, dann $C_{ii}=0$. Andernfalls, wie in der verknüpften Antwort gezeigt, der Eintrag, der enthält$x=0$ ist $-1$. Die resultierende Matrix ist also so

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

Answer 1

Wie der Thread mathematica.se zeigt, ist die Lösung von$$ \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}$$ist nicht einzigartig. Es gibt eine nicht triviale Lösung und die andere Lösung ist$u \equiv 0$.

Formulierung der Gleichung als $$ -\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,$$ wir sehen, dass die Geschwindigkeit der Advektion ist $-\operatorname{sign}(x)$. Das heißt, Masse wird immer in Richtung transportiert$x=0$. Dies erklärt auch die Form der Lösung aus dem mathematica.se-Thread , die bei nicht differenzierbar ist$x=0$.

Nach den üblichen Schritten leiten wir die schwache Form ab $$ \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, $$ was vereinfacht zu $$ -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 $$ unter der Vorraussetzung, dass $u,v$ sind kontinuierlich in $x=0$. Nehmen$u,v \in H^1_0(\Omega)$Dies ist in der Tat aufgrund der Einbettung von Sobolev der Fall.

Wir diskretisieren den Raum $H^1_0(\Omega)$ durch Standardhutfunktionen $\varphi_i$ die auf einem äquidistanten Raster der Größe platziert sind $h$. Das heißt, wir haben$V_h = \operatorname{span}\left\{ \varphi_i : i \in \mathcal{I} \right\} \subset H^1_0(\Omega)$, wo $\mathcal{I}$ ist ein Indexsatz.

Auf dieser Basis konstruieren wir die Matrizen $A$ und $B(w)$, wo $$\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} $$ Hier die Matrix $B$ hängt immer noch von einer Funktion ab $w \in V_h$. Dies führt zu dem (diskreten) Festpunktproblem$$A \vec{u} + B(u_h) \vec{u} = \vec{0},$$ wo $\vec{u}$ bezeichnet die Koordinaten von $u_h \in V_h$.

Wir wenden eine Festpunktiteration an, indem wir das Problem wie folgt linearisieren:

1. Wählen $u_0 \in V_h$ und setzen $n = 0$.
2. Lösen $\displaystyle \left(A + B(u_n)\right) \vec{u}_{n+1} = \vec{0}$ erhalten $\vec{u}_{n+1}$.
3. Überprüfen Sie das Konvergenz- / Stoppkriterium.
4. Wenn das Kriterium nicht erfüllt ist, erhöhen Sie $n$ und fahren Sie mit Schritt 2 fort.

Ich habe dieses Schema im folgenden Python-Skript schnell zusammen gehackt (es ist sehr ineffizient und verwendet nicht einmal spärliche Matrizen). Es konvergiert immer zu$u \equiv 0$, auch wenn sehr nahe an der anderen Lösung begonnen. Eine nicht triviale Lösung kann erhalten werden, wenn eine rechte Seite ungleich Null angewendet wird (auskommentiert).

#!/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()

Ich würde vorschlagen, dass Pseudo-Zeitschritte (oder pseudo-transiente Fortsetzung) von einer anfänglichen Vermutung ungleich Null ausgehen, um die andere nicht triviale Lösung zu berechnen.

Hier ist der Grund (bitte korrigieren Sie mich, wenn ich falsch liege): Wenn wir die Lösung als stationären Zustand der zeitabhängigen Gleichung betrachten, sehen wir, dass der Diffusionsterm (Massenverteilung) den Advektivterm (Transport in Richtung) genau ausgleicht $x=0$). Somit kann im stationären Zustand aufgrund der Randbedingungen und des Strömungsfeldes keine Masse in das System eintreten oder aus diesem austreten. In der Übergangsphase kann die Masse nach Bedarf noch in das System eintreten oder aus dem System austreten, um den stationären Zustand zu erreichen. Daher erscheint mir eine zeitschrittbasierte Methode angemessener als ein Festkomma oder eine Art Newtonsche Iteration.

Für die Festkomma-Iteration vermute ich das $A + B(w)$ ist immer invertierbar, außer für $w \in H^1_0$die nicht triviale Lösung sein. Da wir diese nicht triviale Lösung in nicht genau darstellen können$V_h$, wir enden immer mit $u \equiv 0$. Daher ist eine Festpunktiteration hier nicht geeignet.