Computing the Action of the Resolvent Operator (Stationary PDE)

DrPapa · February 1, 2024, 5:17pm

I am interested in computing the action of the resolvent operator R_A = \left(zI-A\right)^{-1} at a range of z\in\mathbb{C}, where

Ag= F\cdot\nabla g,

and

F:\mathbb{R}^n\rightarrow\mathbb{R}^n

Neumman boundary conditions are needed if F points outward on the boundary. For simplicity I am only considering axis-aligned box domains.

I currently have something working using ApproxFun.jl (Below). However, ApproxFun.jl limits the dimensionality of the problem, so I would like to represent this problem numerically in a form that supports higher dimensions. My thought is to use finite difference for the discretization. However, I am a bit lost on how to proceed w/ forming the discretization.

Looking through the Julia ecosystem I came across DiffEqOperators.jl, SciMLOperators.jl, ParallelStencil.jl, Stencils.jl, and MethodOfLines.jl. Of these, the only one I could figure out how to use was MethodOfLines.jl. However, it produces a NonlinearProblem, which seems unnecessary.

What is the recommended approach to solving such a problem? e.g. sparse matrix, lazy operators, stencils, … Based on other posts, here, my impression is that stencils are the way to go. However, I was surprised that I couldn’t find any packages that had pre-built stencils for such a problem. Also, can stencils be used to solve linear systems?

If you cannot tell already, I know very little about numerically solving PDE’s. Any advice is greatly appreciated. Lastly, if I need to compute the action at different values of z is there a factorization or form better suited for this?

using ApproxFun

F(x) = [2*x[1] - 8*x[1]^3,
           2*x[2] - 8*x[2]^3]
g(x) = 1-x[1]-x[2]^2

domain = Interval(-1.0, 1.0)^2
x̂, ŷ = Fun(domain)
ĝ = g([x̂, ŷ])
F̂ = F([x̂, ŷ])
∂x = Derivative(space(ĝ), [1,0])
∂y = Derivative(space(ĝ), [0,1])

A = F̂[1]*∂x + F̂[2]*∂y # A = F⋅∇
I = IdentityOperator(space(ĝ))

z = 10.0 + 10im
Rg = \(z*I - A, ĝ; tolerance = 1e-3)