# Problem with differential equation in mathematica

I am trying to solve the following differential equation:

$\frac{d^2φ}{dr^2} = -\frac{2}{r}\frac{dφ}{dr} + F(φ,r)$,

where $F(φ,r)$ is a very complicated expression that depends on $φ$ and $r$.

This can't be done analytically, so it must be done numerically.

What I tried to do is the following:

NDSolve[{x''[y] + 2/(y + 0.01)*x'[y] == Y1[x[y]],
x[0] == 0.4,
x'[0] == 0}, x,
{y, 0.1, 100},
Method -> Automatic]


and then I wanted to plot the result. (The 0.01 in the denominator is fine, it's something we use in field theory.)

The problem is that Mathematica gives a solution, which is not what it should be from the theory and I guess it has to do with the following errors that appear when I try to solve the D.E.

Are there any suggestions?

P.S: I should also have mentioned that I am using Mathematica 11.

P.S2: The form of Y1 is

3(-4.23777208123758*^-7 E^(-7.0898154036220635 x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))])
DiracDelta[-4.5 + x] (x Sqrt[1/(1 - (3 Log[x])/(50 \[Pi]))])^(
3/2) (1 + 0.26446386728801075/(
x Sqrt[1/(1 - (3 Log[x])/(50 \[Pi]))])) -
4.23777208123758*^-7 E^(-7.0898154036220635 x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))])
HeavisideTheta[-4.5 + x] (x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))])^(
3/2) (-(0.26446386728801075/(
x^2 Sqrt[1/(1 - (3 Log[x])/(50 \[Pi]))])) - (
0.002525443904885155 Sqrt[1/(1 - (3 Log[x])/(50 \[Pi]))])/
x^2) - 4.23777208123758*^-7 E^(-7.0898154036220635 x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))])
HeavisideTheta[-4.5 + x] (x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))])^(
3/2) (1 + 0.26446386728801075/(
x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))])) (-7.0898154036220635 Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))] -
0.06770275002573074 (1/(1 - (3 Log[x])/(50 \[Pi])))^(3/2)) -
6.35665812185637*^-7 E^(-7.0898154036220635 x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))]) HeavisideTheta[-4.5 + x] Sqrt[
x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))]] (1 + 0.26446386728801075/(
x Sqrt[1/(1 - (3 Log[x])/(50 \[Pi]))])) (Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))] + (
3 (1/(1 - (3 Log[x])/(50 \[Pi])))^(3/2))/(
100 \[Pi])) + ((-HeavisideTheta[-4.5 + x] +
HeavisideTheta[-1.8 + x]) (-0.00003904000000000002 +
0.000035945364096363873 x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))]) (-0.0013364302035827594 Sqrt[
1/(1 - (3 Log[x])/(50 \[Pi]))] -
0.000012761968379850251 (1/(1 - (3 Log[x])/(50 \[Pi])))^(
3/2) + (0.00009600000000000002 x)/(1 - (3 Log[x])/(
50 \[Pi]))^2 + (0.01005309649148734 x)/(
1 - (3 Log[x])/(50 \[Pi]))))/(0.0005219000000000003 -
0.0013364302035827594 x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))] + (0.00502654824574367 x^2)/(
1 - (3 Log[x])/(50 \[Pi])))^2 - ((-DiracDelta[-4.5 + x] +
DiracDelta[-1.8 + x]) (-0.00003904000000000002 +
0.000035945364096363873 x Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))]))/(
0.0005219000000000003 -
0.0013364302035827594 x Sqrt[1/(1 - (3 Log[x])/(50 \[Pi]))] + (
0.00502654824574367 x^2)/(
1 - (3 Log[x])/(
50 \[Pi]))) - ((-HeavisideTheta[-4.5 + x] +
HeavisideTheta[-1.8 + x]) (0.000035945364096363873 Sqrt[1/(
1 - (3 Log[x])/(50 \[Pi]))] +
3.4325294263045497*^-7 (1/(1 - (3 Log[x])/(50 \[Pi])))^(
3/2)))/(0.0005219000000000003 -
0.0013364302035827594 x Sqrt[1/(1 - (3 Log[x])/(50 \[Pi]))] + (
0.00502654824574367 x^2)/(1 - (3 Log[x])/(50 \[Pi]))) +
1/2 (-HeavisideTheta[-1.8 + x] +
HeavisideTheta[
x]) (-0.005671852322897651 x^2 (1/(1 - Log[x]/100))^(3/2) -
0.000028359261614488252 x^2 (1/(1 - Log[x]/100))^(5/2) + (
8.000000000000001*^-6 x)/(1 - (3 Log[x])/(50 \[Pi]))^2 - (
0.000015915494309189537 x^2 ((
0.5026548245743668 x)/(1 - Log[x]/100)^2 + (
100.53096491487337 x)/(1 - Log[x]/100)) (1 - Log[x]/
100))/(1 - (3 Log[x])/(50 \[Pi]))^2 + (
0.0008377580409572782 x)/(1 - (3 Log[x])/(50 \[Pi])) - (
3 x^3 (-5.41 +
Log[(50.26548245743668 x^2)/(1 - Log[x]/100)]))/(
31250 \[Pi] (1 - (3 Log[x])/(50 \[Pi]))^3) - (
2 x^3 (-5.41 +
Log[(50.26548245743668 x^2)/(1 - Log[x]/100)]))/(
625 (1 - (3 Log[x])/(50 \[Pi]))^2)) +
1/2 (-DiracDelta[-1.8 + x] +
DiracDelta[x]) (-0.00002193245422464303 -
0.0018906174409658836 x^3 (1/(1 - Log[x]/100))^(3/2) + (
0.0004188790204786391 x^2)/(1 - (3 Log[x])/(50 \[Pi])) - (
x^4 (-5.41 +
Log[(50.26548245743668 x^2)/(1 - Log[x]/100)]))/(
1250 (1 - (3 Log[x])/(50 \[Pi]))^2))) + (
4 \[Pi] x^4 (27/(25000 \[Pi] x (1 - (3 Log[x])/(50 \[Pi]))^2) + (
39 Sec[1/2 Sqrt[13/
3] (3.67407 + Log[1/(100 (1 - (3 Log[x])/(50 \[Pi])))])]^2)/(
50000 \[Pi] x (1 - (3 Log[x])/(50 \[Pi]))^2) + (
3 Sqrt[39]
Tan[1/2 Sqrt[13/
3] (3.67407 + Log[1/(100 (1 - (3 Log[x])/(50 \[Pi])))])])/(
25000 \[Pi] x (1 - (3 Log[x])/(50 \[Pi]))^2)))/(1 - (3 Log[x])/(
50 \[Pi]))^2 + (
12 x^3 (9/(500 (1 - (3 Log[x])/(50 \[Pi]))) + (
Sqrt[39] Tan[
1/2 Sqrt[13/
3] (3.67407 + Log[1/(100 (1 - (3 Log[x])/(50 \[Pi])))])])/(
500 (1 - (3 Log[x])/(50 \[Pi])))))/(
25 (1 - (3 Log[x])/(50 \[Pi]))^3) + (
16 \[Pi] x^3 (9/(500 (1 - (3 Log[x])/(50 \[Pi]))) + (
Sqrt[39] Tan[
1/2 Sqrt[13/
3] (3.67407 + Log[1/(100 (1 - (3 Log[x])/(50 \[Pi])))])])/(
500 (1 - (3 Log[x])/(50 \[Pi])))))/(1 - (3 Log[x])/(50 \[Pi]))^2

• What is Y1 and there is a syntax problem, why you use brackets for x[y] ? May 24 '17 at 12:42
• Y1 is the F(φ,r). Sorry, I changed notation a bit. If I don't use the brackets, Mathematica says that the function x is with no arguments and does not solve it. May 24 '17 at 13:07
• Your problem seems to be the definition of Y1 therefore the question of @optimalcontrol. When I define Y1[r_]:=r it works like it should without any errors. We need to know Y1 to help you. May 24 '17 at 13:13
• Thank you for your time. I have included the form of Y1 in the main post. May 24 '17 at 13:25
• you will probably make that error go away if you define Y1 to take only numeric arguments. Then you will only throw an error in the unlikely event x hits "exactly" 4.5. (I suspect your usage of delta is ill posed.. ) May 24 '17 at 15:47

As suggested in the comment, you need to define Y1 in a proper way, i.e,

Y1[x_] := paste your long expression here.


But still NDSolve gives an error,

NDSolve::deltad: NDSolve cannot handle discontinuities where the argument of DiracDelta depends on any variable besides the temporal independent variable.

This can tackled by using Method -> {"DiscontinuityProcessing" -> False}.

Once again you will be face with another error,

NDSolve::nlnum:

which show that there are complex numbers, to avoid this you can take Re[Y1[x[y]]],

sol = NDSolve[{x''[y] + 2/(y + 0.01)*x'[y] == Re[Y1[x[y]]], x[0] == 0.4, x'[0] == 0}, x,
{y, 0.1, 100}, Method -> {"DiscontinuityProcessing" -> False}]

Plot[x[y] /. sol, {y, 0.1, 100}]


• notice you get the exact same result if you simplify Y1 by setting all of the DiracDelta 's to zero. May 24 '17 at 16:00
• Thank you for all the help guys. I really appreciate it. I got the same graph as this one, but the thing is that in theory it should not be that way. It should remain zero to infinity and not crossing to negative values. That's the problem. May 24 '17 at 16:21
• Sorry, didn't know about that. And honestly didn't mean to offend you in any way. May 24 '17 at 16:36
• @zhk I did it. I clicked the tick mark and the vote up May 24 '17 at 16:40
• @george2079 Thanks for your nice suggestion.
– zhk
May 25 '17 at 8:34