Elemento finito (1D) per problemi non lineari stazionari

Sep 22 2020

Ho bisogno di risolvere con elementi finiti lineari l'equazione $$\frac{\partial }{\partial x}\Bigl(\text{sgn}(x) u \Big) +\frac{\partial}{\partial x} \Bigl[ \sqrt{u} \frac{\partial u}{\partial x} \Bigr] =0$$

con condizioni al contorno $u(-L)=u(L)=0$ dove $L=6$

(È la versione in stato stazionario dell'equazione qui descritta: diffusione di avvezione non lineare con termine di avvezione non differenziabili )

prendo $v \in H_0^1(-L,L)$ e dopo i soliti passaggi ottengo $$- \int u(x) \text{sgn}(x) v' dx - \int \sqrt{u} u' v' dx = 0$$

Quindi, utilizzando elementi finiti lineari: $$- \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$$

che porta al sistema non lineare (setting$U=[u_0,\ldots,u_N]$)

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

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

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

Ora, voglio risolvere questa equazione non lineare con iterazioni a punto fisso , quindi ho impostato$$CU^{k+1} = -A(U^k)U^k$$ e risolvere iterativamente quei sistemi lineari.

Il problema: sfortunatamente, l'iterazione del punto fisso mi dà NaNe non riesco a trovare la soluzione. È perché il problema è mal posto o ho fatto qualcosa di sbagliato con la mia idea di iterazioni di fixpoint?

Dopo il commento di @ cos_theta, ho modificato il mio codice con la giusta formulazione debole, ma non si riesce ancora a trovare la soluzione. Fondamentalmente, ho creato due funzioni, una in cui assemblo la matrice$A(U)$e l'altro in cui assemblo la matrice $C$. Quindi ho il ciclo di iterazione a virgola fissa.

In particolare, la matrice $A(U)$ corrisponde a $$\int \sqrt{ \sum_k u_k \phi_k } \sum_j u_j \phi_j' \phi_i' dx = 0$$

quindi è tridiagonale e, ad esempio, l'entrata diagonale è $$\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 $$

dove i valori $\sqrt{u_{i-1}}$, $\sqrt{u_i}$, $\sqrt{u_{i+1}}$ sono dati dall'iterazione precedente.

Per la matrice $C$, Ce l'ho $$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$$ Se l'intervallo non contiene $x=0$, poi $C_{ii}=0$. Altrimenti, come mostrato nella risposta collegata, la voce che contiene$x=0$ è $-1$. Quindi la matrice risultante è così

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

Risposte

1 cos_theta Sep 26 2020 at 06:52

Come mostra il thread matematica.se , la soluzione di$$ \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}$$non è unico. C'è una soluzione non banale e l'altra è$u \equiv 0$.

Formulando l'equazione come $$ -\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,$$ vediamo che la velocità dell'avvezione è $-\operatorname{sign}(x)$. Cioè, la massa viene sempre trasportata verso$x=0$. Questo spiega anche la forma della soluzione dal thread matematica.se , che non è differenziabile in$x=0$.

Seguendo i soliti passaggi, deriviamo la forma debole $$ \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, $$ che semplifica a $$ -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 $$ purché $u,v$ sono continui in $x=0$. Prendendo$u,v \in H^1_0(\Omega)$, questo è effettivamente il caso a causa dell'incorporamento di Sobolev.

Discretizziamo lo spazio $H^1_0(\Omega)$ dalle funzioni cappello standard $\varphi_i$ che sono posti su una griglia equidistante di dimensioni $h$. Cioè, abbiamo$V_h = \operatorname{span}\left\{ \varphi_i : i \in \mathcal{I} \right\} \subset H^1_0(\Omega)$, dove $\mathcal{I}$ è un insieme di indici.

Usando questa base, costruiamo le matrici $A$ e $B(w)$, dove $$\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} $$ Qui, la matrice $B$ dipende ancora da qualche funzione $w \in V_h$. Ciò dà origine al problema (discreto) del punto fisso$$A \vec{u} + B(u_h) \vec{u} = \vec{0},$$ dove $\vec{u}$ denota le coordinate di $u_h \in V_h$.

Applichiamo un'iterazione in virgola fissa linearizzando il problema come segue:

Scegliere $u_0 \in V_h$ e impostare $n = 0$.
Risolvere $\displaystyle \left(A + B(u_n)\right) \vec{u}_{n+1} = \vec{0}$ ottenere $\vec{u}_{n+1}$.
Verificare il criterio di convergenza / arresto.
Se il criterio non è soddisfatto, aumentare $n$ e vai al passaggio 2.

Ho rapidamente hackerato questo schema insieme nel seguente script Python (è altamente inefficiente e non usa nemmeno matrici sparse). Converge sempre a$u \equiv 0$, anche se avviato molto vicino all'altra soluzione. Si può ottenere una soluzione non banale se viene applicato un lato destro diverso da zero (commentato).

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

Suggerirei lo pseudo time stepping (o la continuazione pseudo-transitoria) iniziato da un'ipotesi iniziale diversa da zero per calcolare l'altra soluzione non banale.

Ecco perché (correggimi se sbaglio): considerando la soluzione come lo stato stazionario dell'equazione dipendente dal tempo, vediamo che il termine diffusivo (distribuzione della massa) bilancia esattamente il termine advettivo (trasporto verso $x=0$). Pertanto, nello stato stazionario, nessuna massa può entrare o uscire dal sistema a causa delle condizioni al contorno e del campo di flusso. Nella fase transitoria, la massa può ancora entrare o uscire dal sistema secondo necessità per raggiungere lo stato stazionario. Pertanto, un metodo basato sul passaggio del tempo mi sembra più appropriato del punto fisso o di qualche tipo di iterazione di Newton.

Per l'iterazione del punto fisso, lo sospetto $A + B(w)$ è sempre invertibile, ad eccezione di $w \in H^1_0$essendo la soluzione non banale. Poiché non possiamo rappresentare esattamente questa soluzione non banale in$V_h$, finiamo sempre con $u \equiv 0$. Pertanto, l'iterazione in virgola fissa non è adatta qui.