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
$\endgroup$
  • $\begingroup$ What is Y1 and there is a syntax problem, why you use brackets for x[y] ? $\endgroup$ – optimal control May 24 '17 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$ – A_user_with_NoName May 24 '17 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 May 24 '17 at 13:13
  • $\begingroup$ Thank you for your time. I have included the form of Y1 in the main post. $\endgroup$ – A_user_with_NoName May 24 '17 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 May 24 '17 at 15:47
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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.