使用隐式欧拉和共轭梯度线性求解器求解具有非零 Dirichlet BC 的热方程？

Question

许多用户询问如何使用非零 Dirichlet BC 和内部线性求解器的共轭梯度求解热方程 u_t = u_xx。在转向更困难的抛物线 PDE 版本之前，这是一个常见的简化 PDE 问题。 DifferentialEquations.jl 是怎么做到的？

Answer 1

让我们分步解决这个问题。首先，让我们用 Dirichlet BCs 为离散热方程构建线性算子。离散化 can be found on this Wiki page 的讨论表明，中心差分法通过 (u[i-1] - 2u[i] + u[i+1])/dx^2 给出了二阶导数的二阶离散化。这与乘以 [1 -2 1]*(1/dx^2) 的三对角矩阵相同，所以让我们从构建这个矩阵开始：

using LinearAlgebra, OrdinaryDiffEq
x = collect(-π : 2π/511 : π)

## Dirichlet 0 BCs

u0 = @. -(x).^2 + π^2

n = length(x)
A = 1/(2π/511)^2 * Tridiagonal(ones(n-1),-2ones(n),ones(n-1))

注意我们已经隐含地简化了结尾，因为 (u[0] - 2u[1] + u[2])/dx^2 = (- 2u[1] + u[2])/dx^2 当左 BC 为零时，所以该术语从 matmul 中删除。然后我们使用导数的这种离散化来求解热方程：

function f(du,u,A,t)
    mul!(du,A,u)
end

prob = ODEProblem(f,u0,(0.0,10.0),A)
sol = solve(prob,ImplicitEuler())

using Plots
plot(sol[1])
plot!(sol[end])

现在我们让 BCs 非零。请注意，我们只需要加回之前删除的 u[0]/dx^2，所以我们有：

## Dirichlet non-zero BCs
## Note that the operator is no longer linear
## To handle affine BCs, we add the dropped term

u0 = @. (x - 0.5).^2 + 1/12
n = length(x)
A = 1/(2π/511)^2 * Tridiagonal(ones(n-1),-2ones(n),ones(n-1))

function f(du,u,A,t)
    mul!(du,A,u)
    # Now do the affine part of the BCs
    du[1] += 1/(2π/511)^2 * u0[1]
    du[end] += 1/(2π/511)^2 * u0[end]
end

prob = ODEProblem(f,u0,(0.0,10.0),A)
sol = solve(prob,ImplicitEuler())

plot(sol[1])
plot!(sol[end])

现在让我们换掉线性求解器。 The documentation建议你在这里使用LinSolveCG，它看起来像：

sol = solve(prob,ImplicitEuler(linsolve=LinSolveCG()))

这有一些优点，因为它有一个有助于调节的规范处理。但是，文档还指出您可以构建自己的线性求解器例程。这是通过给 Val{:init} 分派 returns 用作线性求解器的类型来完成的，所以我们这样做：

## Create a linear solver for CG
using IterativeSolvers

function linsolve!(::Type{Val{:init}},f,u0;kwargs...)
  function _linsolve!(x,A,b,update_matrix=false;kwargs...)
    cg!(x,A,b)
  end
end

sol = solve(prob,ImplicitEuler(linsolve=linsolve!))

plot(sol[1])
plot!(sol[end])

我们得到了线性求解器的 Krylov 方法（共轭梯度）的非零 Dirichlet 热方程，使其成为 Newton-Krylov 方法。

使用隐式欧拉和共轭梯度线性求解器求解具有非零 Dirichlet BC 的热方程？

Solve the Heat Equation with non-zero Dirichlet BCs with Implicit Euler and Conjugate Gradient Linear Solvers?

julia

differential-equations

differentialequations.jl