1
$\begingroup$

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.

errors

Are there any suggestions?

Thank you very much in advance for your time.

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
Improve this question
$\endgroup$
6
  • $\begingroup$ What is Y1 and there is a syntax problem, why you use brackets for x[y] ? $\endgroup$
    – optimal control
    Commented May 24, 2017 at 12:42
  • $\begingroup$ 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. $\endgroup$
    – user49048
    Commented May 24, 2017 at 13:07
  • 1
    $\begingroup$ 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. $\endgroup$
    – Julien Kluge
    Commented May 24, 2017 at 13:13
  • $\begingroup$ Thank you for your time. I have included the form of Y1 in the main post. $\endgroup$
    – user49048
    Commented May 24, 2017 at 13:25
  • $\begingroup$ 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.. ) $\endgroup$
    – george2079
    Commented May 24, 2017 at 15:47

1 Answer 1

Reset to default
2
$\begingroup$

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}]

enter image description here

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ notice you get the exact same result if you simplify Y1 by setting all of the DiracDelta 's to zero. $\endgroup$
    – george2079
    Commented May 24, 2017 at 16:00
  • $\begingroup$ 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. $\endgroup$
    – user49048
    Commented May 24, 2017 at 16:21
  • $\begingroup$ Sorry, didn't know about that. And honestly didn't mean to offend you in any way. $\endgroup$
    – user49048
    Commented May 24, 2017 at 16:36
  • $\begingroup$ @zhk I did it. I clicked the tick mark and the vote up $\endgroup$
    – user49048
    Commented May 24, 2017 at 16:40
  • $\begingroup$ @george2079 Thanks for your nice suggestion. $\endgroup$
    – zhk
    Commented May 25, 2017 at 8:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.