# GreenFunction Computation for perturbed Laplacian

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!!

• Commented Mar 30, 2021 at 7:06
• Hi Ulrich, Beyond the fact that we both mentioned GreenFunction method, can you tell me how that post has a similar problem to me?
– jcm
Commented Mar 30, 2021 at 20:09
• Hi, it's the same problem , and another view, not more. If you could provide some boundary conditions it might be possible to get numerical approximation. Commented Mar 30, 2021 at 21:04
• Ah, I see! How would a numerical approximation would work if I took {x, y} \[Element] Disk[] with DirichletCondition[u[x, y] == 0, True] ? GreenFunction doesn't seem to like that either
– jcm
Commented Mar 30, 2021 at 22:05

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 ]


• Well, zero Dirichlet b.c. isn't that suitable actually. Just consider the GreenFunction[{-Laplacian[u[x, y], {x, y}]}, u[x, y], {x, y} ∈ FullRegion[2], {m, n}] case, the Green's function isn't tending to zero at infinity. Commented Apr 2, 2021 at 6:45
• @xzcd Thanks, I realize that these finite b.c. doesn't approximate the infinite case, I only tried to show the principle way to get a numerical approximation. Commented Apr 2, 2021 at 10:22

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,

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.

• This approach looks very promising, Thanks! It works perfectly for constant coefficients. For my use case I need it to work for non-constant coefficients. For instance, a(x,y) = x^2 instead of simply a constant 2. Looks like DChange has trouble with this, as I am getting an 'iteration limit exceeded' . Do you have much experience with that sort of error?
– jcm
Commented Apr 1, 2021 at 20:51
• @jcm After a second thought I notice the DChange-based solution is wrong. (I forgot the RHS of the equation. ) Answer rewritten, now it should be correct. Sadly I can't think out a way to generalize this method to handle the new example you mentioned. Commented Apr 5, 2021 at 10:10