非线性 Poisson 方程求解

非线性方程在实际问题中经常出现,这里详细介绍求解非线性 Poisson 方程的典型方法,如

  • Newton-Galerkin 方法
    • 先应用 Newton 方法线性化连续的弱形式,再用有限元离散次迭代求解
    • 先应用有限元离散连续的弱形式,再应用 Newton 方法迭代求解
  • Picard 迭代方法
    • 又称定点迭代

Newton-Galerkin 方法

首先给出一个扩散系数为非线性的例子. 求解区域记为 , 边界 .

满足如下的边界条件:

在 Poisson 方程两端分别乘以测试函数 , 利用分部积分,可得到其对应的连续弱形式

的一个逼近,记 , 代入连续弱形式

其中 处 Taylor 展开,可得

替换连续弱形式中的 , 并忽略掉其中 可得

给定求解区域 上的网格离散 , 构造 维的有限维空间 , 其 $N$ 个全局基函数组成的行向量函数记为

对于有限元程序设计实现来说,并不会显式构造出全局基函数,实际基函数的求值计 算都发生网格单元或网格单元的边界上。设每个网格单元 局部基函数个数为 个,其组成的行向量函数记为

所有基函数梯度对应的向量

其中

对应的单元矩阵为

对应的单元矩阵为

对应的单元列向量为

下面讨论边界条件相关的矩阵和向量组装。设网格边界边(2D)或边界面(3D)上的局部基函数个数为 $m$ 个,其组成的行向量函数记为

是一个边界边或边界面,则 对应的矩阵为

对应的向量为

对应的向量为

基于 FEALPy 的程序实现

设求解区域为 真解设为

非线性扩散系数设为

边界条件设为纯 Dirichlet 边界条件,其线性化的连续弱形式为

首先导入必要的模块

import numpy as np

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve

from fealpy.decorator import cartesian, barycentric
from fealpy.mesh import MeshFactory as MF
from fealpy.functionspace import LagrangeFiniteElementSpace
from fealpy.boundarycondition import DirichletBC
from fealpy.tools.show import showmultirate

接着定义模型数据函数

@cartesian
def solution(p):
    # 真解函数
    pi = np.pi
    x = p[..., 0]
    y = p[..., 1]
    return np.cos(pi*x)*np.cos(pi*y)

@cartesian
def gradient(p):
    # 真解函数的梯度
    x = p[..., 0]
    y = p[..., 1]
    pi = np.pi
    val = np.zeros(p.shape, dtype=np.float64)
    val[...,0] = -pi*np.sin(pi*x)*np.cos(pi*y)
    val[...,1] = -pi*np.cos(pi*x)*np.sin(pi*y)
    return val #val.shape ==p.shape

@cartesian
def source(p):
    # 源项
    x = p[...,0]
    y = p[...,1]
    pi = np.pi
    val = 2*pi**2*(3*np.cos(pi*x)**2*np.cos(pi*y)**2-np.cos(pi*x)**2-np.cos(pi*y)**2+1)*np.cos(pi*x)*np.cos(pi*y)
    return val

@cartesian
def dirichlet(p):
    return solution(p)

实现 $(a_u'(u^0)\nabla u^0\cdot\delta u, \nabla v)$ 对应的组装代码

def nolinear_matrix(uh, q=3):

    space = uh.space
    mesh = space.mesh

    qf = mesh.integrator(q, 'cell')
    bcs, ws = qf.get_quadrature_points_and_weights()

    cellmeasure = mesh.entity_measure('cell')

    cval = 2*uh(bcs)[...,None]*uh.grad_value(bcs) # (NQ, NC, GD)
    phi = space.basis(bcs)       # (NQ, 1, ldof)
    gphi = space.grad_basis(bcs) # (NQ, NC, ldof, GD)

    B = np.einsum('q, qcid, qcd, qcj, c->cij', ws, gphi, cval, phi, cellmeasure)

    gdof = space.number_of_global_dofs()
    cell2dof = space.cell_to_dof()
    I = np.broadcast_to(cell2dof[:, :, None], shape=B.shape)
    J = np.broadcast_to(cell2dof[:, None, :], shape=B.shape)
    B = csr_matrix((B.flat,(I.flat,J.flat)), shape=(gdof,gdof))

    return B

准备好求解参数

p = 1 # 有限元空间次数
tol = 1e-8 # 非线性迭代误差限

构造 $\Omega=[0, 1]^2$ 上的初始三角形网格,$x$ 和 $y$ 方向都剖分 10 段

domain = [0, 1, 0, 1]
mesh = MF.boxmesh2d(domain, nx=10, ny=10, meshtype='tri')

构造 mesh 上的 $p$ 次 Lagrange 有限元空间 space, 并定义一个该空间的函数 u0du, 其所有自由度系数默认取 0

space = LagrangeFiniteElementSpace(mesh, p=p) # p 的线性元,
u0 = space.function() # 有限元函数 u0 = 0
du = space.function() # 有限元函数 du = 0

设置 Dirichlet 边界条件

isDDof = space.set_dirichlet_bc(dirichlet, u0)
isIDof = ~isDDof

计算载荷向量

b = space.source_vector(source)

定义扩散系数函数

@barycentric
def dcoefficient(bcs):
    return 1+u0(bcs)**2

非线性迭代求解

while True:
    A = space.stiff_matrix(c=dcoefficient)
    B = nolinear_matrix(u0)
    U = A + B
    F = b - A@u0
    du[isIDof] = spsolve(U[isIDof, :][:, isIDof], F[isIDof]).reshape(-1)
    u0 += du
    err = np.max(np.abs(du))
    print(err)
    if err < tol:
       break

第二个是反应项为非线性的例子

满足如下的边界条件:

results matching ""

    No results matching ""