许多用户询问如何使用非零 Dirichlet BC 和内部线性求解器的共轭梯度求解热方程,u_t = u_xx。在转向更困难的抛物线 PDE 版本之前,这是一个常见的简化 PDE 问题。这是如何在 DifferentialEquations.jl 中完成的?
最佳答案
让我们分步解决这个问题。首先,让我们用 Dirichlet 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])
![enter image description here](https://i.sstatic.net/CQfi3.png)
现在我们使 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])
![enter image description here](https://i.sstatic.net/RpTAa.png)
现在让我们换掉线性求解器。 The documentation建议您应该使用
LinSolveCG
在这里,看起来像:sol = solve(prob,ImplicitEuler(linsolve=LinSolveCG()))
这有一些优点,因为它具有有助于调节的规范处理。但是,文档还指出您可以构建自己的线性求解器例程。这是通过给出 Val{:init}
来完成的。 dispatch 返回用作线性求解器的类型,所以我们这样做:## 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 方法(共轭梯度),使其成为 Newton-Krylov 方法。
关于julia - 用隐式欧拉和共轭梯度线性求解器求解非零狄利克雷 BC 的热方程?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/54545738/