# Could anyone help me to solve the following differential equation?

I have a problem with the solution of this PDE, how can I start to solve it?

$\left(\partial_x^2 + \frac{1}{x}\partial_x + \partial_y^2 \right) F(x,y) = \frac{4\pi I}{x} \left(\frac{\sqrt{x^2+(y\pm M)^2}-2M} {\sqrt{x^2+(y\pm M)^2}+2M}\right)^{1/2}\delta(x-x_0)\delta(y-y_0)$

where $I$ and $M$ are the positive constants.

• unreadable and uncopiable to Mathematica! please copy and paste your mathematica code you rote to solve the equation. when you paste, each line of your code must start with 4 spaces. – Navid Rajil Feb 5 '18 at 6:13
• Are you trying to solve this with the Mathematica software, and if so what have you already tried? If you are only interested in the solution, you may wish to post it on math.stackexchange.com since this site is specifically for questions regarding the Mathematica software. It is rather unlikely that this PDE has a nice closed form solution regardless, so you may have better luck defining boundary conditions and using numerical methods. – eyorble Feb 5 '18 at 7:18
• Welcome to Mathematica.SE. Are you sure you are posting on the right site? There is nothing in your question making it clear that it is concerned with Mathematica software. – m_goldberg Feb 5 '18 at 7:47
• I think that he needs some directions on how to solve this D.E using Mathematica. This is what I understand at least. – A_user_with_NoName Feb 5 '18 at 9:04

This problem is probably more mathematical physics than Mathematica, but it is too interesting for me to ignore. It would help to have some boundary conditions specified in the question, but I will choose some that make the problem doable. Since we are dealing with Cartesian coords, consider a rectangular boundary with 0 <=x<=a and 0 <=y<=b, with D[F[x,y],x] = 0 at x = 0 and F = 0 on the other boundaries. We can see that the RHS of the pde is 0 everywhere except at {x,y} == {x0,y0}, so we can replace x and y under the radical with x0 and y0, which are within the rectangular boundary. Separate variables for the homogeneous eq only to get the form of the solution.

pde = D[F[x, y], x, x] + (1/x)*D[F[x, y], x] + D[F[x, y], y, y] ==
((4*Pi*II)/x)*Sqrt[(Sqrt[x^2 + (y + M)^2] - 2*M)/(Sqrt[x^2 + (y + M)^2] + 2*M)]*DiracDelta[x - x0]*
DiracDelta[y - y0];

F[x_, y_] = X[x] Y[y];

pde[]/F[x, y] == 0 // Expand
(*Derivative[X][x]/X[x] + Derivative[X][x]/(x*X[x]) + Derivative[2[Y][y]/Y[y] == 0*)


The x part and y part must each be separately equal to a constant that could be positive, negative or zero. I choose

xeq = Derivative[X][x]/X[x] + Derivative[X][x]/(x*X[x]) == -alpha^2

DSolve[xeq, X[x], x]
{{X[x] -> C BesselJ[0, alpha x] + C BesselY[0, alpha x]}}


We won't consider the BesselY since it is not well behaved at x = 0. We don't really use the above as the solution, but we will use its form to get a solution. First, expand the DiracDelta[x-x0]/x on the RHS of the pde as

xddeq = (1/x)*DiracDelta[x - x0] == Sum[An*BesselJ[0, (xn*x)/a], {n, 1, Infinity}];


Where xn is the nth root of BesselJ[0,x] = 0. This ensures F = 0 at x = a. We get a warning about no convergence from MMA, but that is the nature of a DiracDelta fn. From orthogonality solve for An by multiplying each side by x BesselJ[0,(xn x)/a] and integrating from 0 to a.

\$Assumptions = a > 0 && b > 0 && 0 < x0 < a && 0 < y0 < b && xn > 0

Integrate[(1/x)*DiracDelta[x - x0]*x*BesselJ[0, (xn*x)/a], {x, 0, a}] ==
An*Integrate[x*BesselJ[0, (xn*x)/a]^2, {x, 0, a}]

Solve[%, An] // Flatten;

An = An /. % /. BesselJ[0, xn] -> 0
(*(2*BesselJ[0, (x0*xn)/a])/(a^2*BesselJ[1, xn]^2)*)


We used the fact that xn is a root of J0. Now

xddeq
(*DiracDelta[x - x0]/x == Sum[(2*BesselJ[0, (x*xn)/a]*BesselJ[0,(x0*xn)/a])/
(a^2*BesselJ[1, xn]^2), {n, 1, Infinity}]*)


Now expand F as

F == (2*Sum[((BesselJ[0, (x*xn)/a]*BesselJ[0, (x0*xn)/a])*f[y])/BesselJ[1, xn]^2, {n, 1, Infinity}])/a^2;


Leave out the summation for now to make simplification easier:

F[x_, y_] = (2/a^2)*((BesselJ[0, (x*xn)/a]*BesselJ[0, (x0*xn)/a])/BesselJ[1, xn]^2)*f[y]

(pde[] // FullSimplify) == pde[]

% /. BesselJ[0, (x xn)/a] BesselJ[0, (x0 xn)/a] -> a^2/2 BesselJ[1, xn]^2 DiracDelta[x - x0]/x

%/(DiracDelta[x - x0]/x) // Apart
(*Derivative[f][y] - (xn^2*f[y])/a^2 == 4*Pi*II*DiracDelta[y - y0]*
Sqrt[(Sqrt[(M + y)^2 + x^2] - 2*M)/(Sqrt[(M + y)^2 + x^2] + 2*M)]*)


The y equation becomes

fyeq = Derivative[f][y] - (xn^2*f[y])/a^2 == 4*Pi*II*DiracDelta[y - y0]*
(Sqrt[(Sqrt[(M + y)^2 + x^2] - 2*M)/(Sqrt[(M + y)^2 + x^2] + 2*M)] /. {x -> x0, y -> y0})


Look at the homogeneous part

DSolve[fyeq[] == 0, f[y], y] // Flatten


We get exponentials, but for our bc's we need trigs

(%[[1, 2]] // ExpToTrig // Simplify) /. {C - C -> A,
C + C -> B}
(*A*Sinh[(xn*y)/a] + B*Cosh[(xn*y)/a]*)


Use the above to form one function for y < y0, and another for y > y0. We have to throw the Cosh term out for the fn that includes y = 0.

f1[y_] = Sinh[(xn y)/a];
f2[y_] = A Sinh[(xn y)/a] + B Cosh[(xn y)/a];


At y = 0, f1==0 is already satisfied, but for f2 we need f2[b] == 0. Solve for A in terms of B

A = A /. Solve[f2[b] == 0, A][];

f2[y_] = f2[y] // Simplify
(*B*Csch[(b*xn)/a]*Sinh[(xn*(b - y))/a]*)


We will form f[y] = f1[ylesser]*f2[ygreater] where ylesser is the lesser of y and y0 and ygreater is the greater of y and y0. The discontinuity in slope at y = y0 equals D[f[y],y]outer - D[f[y],y]inner which turns out to be the Wronskian of the inner and outer functions. So we have.

f1[y] D[f2[y], y] - f2[y] D[f1[y], y] == fyeq[]/DiracDelta[y - y0] // Simplify

B = B /. Solve[%, B][]


and

f[y_] = f1[ylesser] f2[ygreater] // Simplify
(*-((4*Pi*a*II*Sqrt[(Sqrt[(M + y)^2 + x^2] - 2*M)/(Sqrt[(M + y)^2 + x^2] + 2*M)]*
Csch[(b*xn)/a]*Sinh[(xn*ylesser)/a]*Sinh[(xn*(b - ygreater))/a])/xn)*)


Put everything together adding back the summation we get

F[x_, y_] = -4*Pi*a*II*Sqrt[(Sqrt[(M + y0)^2 + x0^2] - 2*M)/(Sqrt[(M + y0)^2 + x0^2] + 2*M)]*(2/a^2)*
Sum[((BesselJ[0, (x*xn)/a]*BesselJ[0, (x0*xn)/a])/BesselJ[1, xn]^2)*
((Csch[(b*xn)/a]*Sinh[(xn*ylesser)/a]*Sinh[(xn*(b - ygreater))/a])/xn), {n, 1, Infinity}]


which is the solution. Plug in some numbers for plotting.

a = 1;
b = 2;
II = 1;
M = 1;
x0 = a/5;
y0 = 2 b/3;


Get the terms that include ylesser and ygreater into a form that we can plot.

g[y_] = Piecewise[{{Sinh[(xn[n] y)/a] Sinh[(xn[n] (b - y0))/a],
y <= y0}, {Sinh[(xn[n] y0)/a] Sinh[(xn[n] (b - y))/a], y > y0}}]


and xn are the roots of J0

xn[n_] = N[BesselJZero[0, n]];


For plotting we can define F this way

F[x_, y_, m_] := (-(1/a))*8*Pi*II*Sqrt[(Sqrt[(M + y0)^2 + x0^2] - 2*M)/(Sqrt[(M + y0)^2 + x0^2] + 2*M)]*
Sum[(Csch[(b*xn[n])/a]*BesselJ[0, (x*xn[n])/a]*BesselJ[0, (x0*xn[n])/a]*g[y])/(xn[n]*BesselJ[1, xn[n]]^2),
{n, 1, m}]

Plot3D[Evaluate[F[x, y, 100]], {x, 0, a}, {y, 0, b},
AxesLabel -> {"x", "y", "F"}, PlotRange -> All] and some of the more interesting cross sections

Presumably it would take many more terms to get an accurate value for the discontinuity peak.