Решить уравнение теплопроводности с ненулевыми БК Дирихле с неявным линейным решателем Эйлера и сопряженного градиента? - PullRequest
Новая
       83

Решить уравнение теплопроводности с ненулевыми БК Дирихле с неявным линейным решателем Эйлера и сопряженного градиента?

0 голосов
/ 06 февраля 2019

Многие пользователи спрашивают, как решить уравнение теплопроводности, u_t = u_xx, с ненулевыми BC Дирихле и с сопряженными градиентами для внутреннего линейного решателя.Это общая упрощенная проблема PDE перед переходом к более сложным версиям параболических PDE.Как это сделать в DifferentialEquations.jl?

1 Ответ

0 голосов
/ 06 февраля 2019

Давайте решим эту проблему поэтапно.Во-первых, давайте построим линейный оператор для дискретного уравнения теплопроводности с БК Дирихле.Обсуждение дискретизации можно найти на этой вики-странице , которая показывает, что метод центральной разности дает дискретизацию второго порядка второй производной на (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 равенноль, поэтому термин исключен из матмул.Затем мы используем эту дискретизацию производной для решения уравнения теплопроводности: 

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

Теперь мы делаем БК ненулевыми.Обратите внимание, что нам просто нужно добавить обратно 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

Теперь давайте поменяем местамилинейный решатель. В документации предлагается, чтобы вы использовали здесь LinSolveCG, который выглядит следующим образом: 

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

В этом есть некоторые преимущества, так как он имеет обработку норм, которая помогает формированию.Однако в документации также говорится, что вы можете создать свою собственную программу линейного решателя.Это делается путем предоставления 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])

И вот мы, ненулевое уравнение теплопроводности Дирихле с методом Крылова (сопряженные градиенты) для линейного решателя, что делает его методом Ньютона-Крылова.

...