FEniCS教程：非线性与时间相关问题求解

# FEniCS 教程：非线性与时间相关问题求解 ## 1. 非线性变分问题的求解 ### 1.1 牛顿法在 PDE 层面的应用 在求解非线性偏微分方程（PDE）时，可将含未知量 $\delta u$ 的项移到等式左侧，得到： $\nabla\cdot\left\{q(u_k)\nabla\delta u\right\}+\nabla\cdot\left\{q^{\prime}(u_k)\delta u\nabla u_k\right\}=-\nabla\cdot\left\{q(u_k)\nabla u_k\right\}$ 将该 PDE 乘以测试函数 $v$ 并在区域 $\Omega$ 上积分，同时对二阶导数进行分部积分，得到其弱形式： $\int_{\Omega}\left\{q(u_k)\nabla\delta u\cdot\nabla v + q^{\prime}(u_k)\delta u\nabla u_k\cdot\nabla v\right\}dx = -\int_{\Omega}q(u_k)\nabla u_k\cdot\nabla v dx$ 变分问题可表述为：找到 $\delta u\in V$，使得对于所有 $v\in\hat{V}$，有 $a(\delta u, v) = L(v)$，其中： $a(\delta u, v) = \int_{\Omega}\left\{q(u_k)\nabla\delta u\cdot\nabla v + q^{\prime}(u_k)\delta u\nabla u_k\cdot\nabla v\right\}dx$ $L(v) = -\int_{\Omega}q(u_k)\nabla u_k\cdot\nabla v dx$ 函数空间 $V$ 和 $\hat{V}$ 与线性泊松问题中的相同。我们需要提供一个初始猜测值，例如求解 $q(u) = 1$ 时的 PDE 解。对应的弱形式 $a_0(u_0, v) = L_0(v)$ 为： $a_0(u, v) = \int_{\Omega}\nabla u\cdot\nabla v dx$ $L(v) = 0$ 之后，进入循环求解 $a(\delta u, v) = L(v)$ 得到 $\delta u$，并计算新的近似值 $u_{k + 1} = u_k + \delta u$。需要注意的是，$\delta u$ 是一个修正项，如果 $u_0$ 在边界的某部分 $\Gamma_D$ 上满足给定的狄利克雷条件，则必须要求 $\delta u$ 在 $\Gamma_D$ 上为 0。 ### 1.2 直接求解非线性变分问题 FEniCS 中的工具可以自动完成之前手动计算和实现的皮卡或牛顿方法。可以使用以下代码直接定义和求解非线性变分问题： ```python problem = NonlinearVariationalProblem(F, u, bcs, J) solver = NonlinearVariationalSolver(problem) solver.solve() ``` 其中，$F$ 对应非线性形式 $F(u; v)$，$u$ 是未知的函数对象，$bcs$ 表示本质边界条件（狄利克雷边界条件对象列表），$J$ 是 $F$ 的雅可比矩阵的变分形式。 #### 1.2.1 定义 $F$ 形式 假设 $q(u)$ 被编码为 Python 函数，$F$ 形式可以通过以下两种方式定义： 方式一： ```python v = TestFunction(V) u_ = Function(V) # 未知量 F = inner(q(u_)*nabla_grad(u_), nabla_grad(v))*dx ``` 方式二： ```python u = TrialFunction(V) v = TestFunction(V) F = inner(q(u)*nabla_grad(u), nabla_grad(v))*dx u_ = Function(V) # 最近计算的解 F = action(F, u_) ``` #### 1.2.2 计算雅可比矩阵 $J$ 雅可比矩阵 $J$ 是 $F$ 在 $u = u^-$ 处沿 $\delta u$ 方向的 Gateaux 导数。可以通过手动计算或使用 UFL 语言的 `derivative` 函数自动计算。 手动计算时，以 $du$ 为试验函数，$J$ 的定义如下： ```python du = TrialFunction(V) v = TestFunction(V) u_ = Function(V) # 最近计算的解 F = inner(q(u_)*nabla_grad(u_), nabla_grad(v))*dx J = inner(q(u_)*nabla_grad(du), nabla_grad(v))*dx + \ inner(Dq(u_)*du*nabla_grad(u_), nabla_grad(v))*dx ``` 使用 `derivative` 函数自动计算时： ```python du = TrialFunction(V) v = TestFunction(V) u_ = Function(V) # 最近计算的解 F = inner(q(u_)*nabla_grad(u_), nabla_grad(v))*dx J = derivative(F, u_, du) # Gateaux 导数，沿 du 方向 ``` #### 1.2.3 设置牛顿法参数 以下代码展示了如何定义非线性变分问题和相关求解器，并设置牛顿法的关键参数： ```python problem = NonlinearVariationalProblem(F, u_, bcs, J) solver = NonlinearVariationalSolver(problem) prm = so ```