解决散热式与隐式Euler和共轭梯度线性解法非零狄氏的BC？

Question

许多用户询问如何解决热方程，u_t = u_xx，具有非零狄利克雷BC和用共轭梯度为内部线性解算器。这是移动到抛物线PDE的更困难的版本之前的共同简化PDE问题。这是如何在DifferentialEquations.jl做了什么？

Answer 1

让我们来解决这个问题的步骤。首先，让我们建立了离散热方程的线性算子与狄氏的BC。离散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])

现在我们做的BC非零。请注意，我们只需要添加回到我们先前下车的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()))

有一些优势，这一点，因为它有一个规范的处理，有助于调理。 Howerver，该文档还指出，你可以建立自己的线性解算器例程。这是通过给予Val{:init}调度返回的线性解算器所用的类型做了，所以我们做的：

## 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法（共轭梯度）为线性求解，使其成为一个牛顿的Krylov法。