4

In 1D the convolution of the unit box function and the modified Bessel function of the second kind $K_0(x)$ works very well.

Clear[f, g, h];
f[x_] := UnitBox[x];
g[x_] := BesselK[0, Abs[x]];
h = Convolve[f[y], g[y], y, x]
Plot[{f[x], g[x], h}, {x, -2, 2}, PlotLegends -> "Expressions"]

enter image description here

However, when I try the same in 2D, the convolution doesn't work.

Clear[f, g, h];
f[x1_, x2_] := UnitBox[x1, x2];
g[x1_, x2_] := BesselK[0, Sqrt[x1^2 + x2^2]];
h = Convolve[f[y1, y2], g[y1, y2], {y1, y2}, {x1, x2}]

How can I solve this problem in 2D? I am only interested in the analytical solution.

Pavlo Fesenko
  • 368
  • 1
  • 10

2 Answers2

1

When you say the 1D convolution of the unit box function and the modified Bessel function of the second kind $K_0(x)$ works very well, I'm not sure what you mean—nor the intended application as $K_0(x)$ is complex for $x<0$.

Note you get a different answer if you restrict the range to [-1,1] in the explicit integral

Integrate[BesselK[0, x - y], {y, -1, 1}]

And this answer is just what you get from the indefinite integral.

Finally, a natural alternative to the 2D unit box is a circle (which generalises the unit interval) , in which case you could use rotational symmetry to obtain the answer.

TheDoctor
  • 2,832
  • 1
  • 14
  • 17
  • I am trying to solve the 2D screened Poisson equation for the unit box source function $f(x_1, x_2)$. The modified Bessel function of the second kind $K_0(|x|)$ is the fundamental solution of this equation. I need to convolve it with the source function in order to get the actual solution. By "works well", I meant that Convolve[] gives me an output which I can use for calculations and plotting (the green line on the plot). My actual problem requires the 2D unit box so a circle is not desirable. I agree, however, that it can be used as an approximation. – Pavlo Fesenko Mar 17 '16 at 15:43
  • I have corrected my question by writing BesselK[0, Abs[x]] in 1D instead of just BesselK[0, x]. So there should be no problems with $x < 0$. In 2D it has already been taken into account BesselK[0, Sqrt[x1^2 + x2^2]]. I am still wondering why the same convolution cannot be calculated in 2D? – Pavlo Fesenko Mar 18 '16 at 09:52
  • Have you considered the possibility that Mathematica just doesn't know a closed form, @Pavlo? In which case, can you not use a numerical method instead? – J. M.'s missing motivation Mar 18 '16 at 10:19
  • @J.M. I agree that it could very likely be the case. On the other hand, is there a workaround to obtain the analytical solution? I have already solved this problem numerically but it's quite difficult to analyze the influence of different parameters on the final result. – Pavlo Fesenko Mar 18 '16 at 12:54
1

Perhaps this suffices:

f[x_, y_] := 
 NIntegrate[
  UnitBox[s, t] BesselK[0, 
    Sqrt[(x - s)^2 + (y - s)^2]], {s, -Infinity, 
   Infinity}, {t, -Infinity, Infinity}]

Visualizing:

Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotRange -> Full, 
 Mesh -> False, PlotPoints -> 25]

Plot could be improved but I do not have time at present.

enter image description here

ubpdqn
  • 60,617
  • 3
  • 59
  • 148
  • Nice idea to use a numerical method to get an idea what the solution looks like. But shouldn't it be y-t in the Sqrt expression inside NIntegrate? Also visually i would expect the solution to resemble more a square than a diagonal line. – Thies Heidecke Jan 08 '18 at 00:11