GreenFunction Computation for perturbed Laplacian

Question

I am happy computing the Green function for the Laplacian

Gsol := GreenFunction[{-Laplacian[u[x, y], {x, y}]}, 
  u[x, y], {x, y} ∈ FullRegion[2], {m, n}]

it gives an analytical solution - which is great.

However, if I slightly perturb my operator

Gsol := GreenFunction[{D[u[x, y], {x, 2}] + 2 D[u[x, y], {y, 2}]}, 
  u[x, y], {x, y} ∈ FullRegion[2], {m, n}]

where I have simply put a 2 infront of the second derivative.

Mathematica is no longer able to provide an analytical or even a numerical approximation to the Green function for this operator. If GreenFunction is not able to handle such an operator, is there another method I can use to do this?

Thanks!!

Answer 1

Here my idea (see comments) to get a numerical approximation for Green's function:

First one needs an approximation of DiracDelta, for example

eps = .001; (*sufficient small*)
dirac = Function[x, Exp[-(x^2/(2 eps))]/Sqrt[2 Pi eps]]      

GreenFunction is the solution of D[u[x, y], {x, 2}] + 2 D[u[x, y], {y, 2}] == dirac[x - \[Xi]] dirac[y - \[Eta]]

 green = ParametricNDSolveValue[{D[u[x, y], {x, 2}] + 2 D[u[x, y], {y, 2}] == dirac[x - \[Xi]] dirac[y - \[Eta]], 
 DirichletCondition[u[x, y] == 0, True]},u , {x, y} \[Element] Disk[] , {\[Xi], \[Eta]}]   

green's function g[\[Xi], \[Eta]][x,y]

Plot3D[g[0.1, 0.3][x, y], {x, y} \[Element] Disk[], PlotRange -> All ]

Answer 2

As a new introduced function, GreenFunction is still fragile in my view. (The post linked by Ulrich in the comment is another example. ) As to your specific problem, the only work-around I can think out at the moment is to directly solve the PDE with Fourier transform, I'll use ft to facilitate coding:

(* Definition of ft isn't included in this post, 
   please find it in the link above. *)

eq = D[u[x, y], {x, 2}] + 2 D[u[x, y], {y, 2}] == - DiracDelta[x - m] DiracDelta[y - n];

rule = HoldPattern@FourierTransform[a_, __] :> a;

teq = ft[eq, x, w1] /. rule /. u -> (U@#2 &)

tteq = ft[teq, y, w2] /. rule /. U[y] -> υ

ttsol = Solve[tteq, υ] // Values // Flatten // First

sol = InverseFourierTransform[ttsol, {w2, w1}, {y, x}]

Notice u -> (U@#2 &) and U[y] -> υ aren't actually necessary, they're just to make the code more readable.

The sol is a bit lengthy and not that easy to simplify, so I'd like to show another way to calculate the inverse Fourier transform that produces a simpler result. We know inverse Fourier transform is essentially an integral, but in this case the generic function has involved in so we can't use Integrate to calculate the inverse transform directly, nevertheless, utilizing the differentiation property of Fourier transform:

FourierTransform[f'[t], t, ω]
(* -I ω FourierTransform[f[t], t, ω] *)

we can calculate derivative of the fundamental solution:

dsol = -(1/Sqrt[2 Pi])^2 Integrate[
   Exp[-I ( w1 x + w2 y)] ttsol w1 w2, {w1, -Infinity, Infinity}, {w2, -Infinity, 
    Infinity}, Assumptions -> (x | y | m | n) ∈ Reals]

(* ConditionalExpression[(
 Sqrt[2] (m - x) Abs[n - y] Sign[n - y])/(π (2 (m - x)^2 + (n - y)^2)^2), n != y] *)

Finally integrate with respect to $x$ annd $y$:

solalter = 
 Integrate[dsol, x, y, Assumptions -> (x | y | m | n) ∈ Reals && n != y] // 
  FullSimplify
(* -(Log[2 (m - x)^2 + (n - y)^2]/(4 Sqrt[2] π)) *)

"Wait, haven't you ignored the constant? " Indeed, but it doesn't hurt because fundamental solution is not unique. Actually it's not hard to show by numeric tests that:

sol - solalter == (-2 EulerGamma + Log[2])/(4 Sqrt[2] π)

Sadly, this method doesn't seem to be easy to generalize to cases like D[u[x, y], {x, 2}] + x^2 D[u[x, y], {y, 2}] == - DiracDelta[x - m] DiracDelta[y - n] as asked in comment.