Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

This is my first question on this site, I hope it is written in an understandable way.

Let's start. My aim is to plot a certain parametric region, with $x$ and $y$ coordinates given by the following formulas. First of all I write down some functions that, once I define a particular potential $U$ of two variables #1 and #2, will calculate some of its derivatives and things like that (since I will need them later for the construction of the plot). Here below are the formulas, the last three are the ones that will be appearing in the following lines of code:

epsilontemp[U_] = D[U[#1,#2],#1])^2/(U[#1,#2])^2/2&;

etatemp[U_] = D[U[#1,#2],{#1,#2}])/U[#1,#2]&;

inttemp[U_] = D[U[#1,#2],#1])/(U[#1,#2])&;

epsilonhelp[U_,a_] = epsilontemp[U][a,1];

etahelp[U_,a_] = etatemp[U][a,1];

inthelp[U_,a_] = inttemp[U][a,1];

Now I write down the potential: it is an infinite series, but Mathematica recognizes it as the sum of two PolyLog functions:


in the following equations I will call always #1 as $a$ and #2 as $b$.

Now: I must find the value of $\frac{a}{b}$ such that epsilontemp[U3][a,b]=1 (this equation can be written more simply in the way I did in the code right below), and I found that this piece of code does the job perfectly (the particular starting point was chosen after I did a "Plot of the equation", in order to find approximately the root):

phiendtemp[b_?NumericQ] := FindRoot[inthelp[U3,x]-Sqrt[2]*b,{x,48/10},MaxIterations->Infinity,Compiled->False][[1,2]];

and in order to work more easily with this quantity, I defined an InterpolatingFunction in the following way:

phiend[b_] = FunctionInterpolation[phiendtemp[btemp],{btemp,3,5}][b];

Now the next step consist in solving numerically a differential equation, with two parameters: one enters simply in the equation while the other enters in the initial condition. The new command of Mathematica ParametricNDSolve was used:

phicmbtemp1 = ParametricNDSolve[y'[eFolds]-b^{-2}*inthelp[U3,y[eFolds]]==0,y[0]==j,y,{eFolds,30,70},{b,j},Compiled->False]

phicmbtemp2[b_,j_,eFolds_] = y[b,j][eFolds]/.phicmbtemp1;

The last step is use the function I defined previously as the initial condition for the differential equation, and I did that in the following way.

phicmb[b_,eFolds_] = phicmbtemp2[b,phiend[b],eFolds];

Up to now the code works perfectly, I also tried to plot phicmb[b,eFolds] for $b=4$ and for $eFolds$ from 30 to 70 (the values I want) and it gives me the results that I expect from the "non numerical" analysis of the equations.

Now I turn to the ParametricPlot: the functions that will be the $x$ and $y$ coordinates are the following (the parameters of the region will be $b$ and $eFolds$):

nsplot[U_,b_,eFolds_] = 1 + 2*b^{-2}*etahelp[U,phicmb[b,eFolds]] - 6*b^{-2}*epsilonhelp[U,phicmb[b,eFolds]];

After I wrote this, I can at last explain what's my problem. The problem is that if I try to calculate things like (in the actual plots, I will use $c$ for the quantity $b$, and $n$ for the quantity $eFolds$)


Mathematica returns

0.85485 + 0.i

that (apart the imaginary part, that is a problem on its own) is compatible with what I expect from a "non numerical" analysis of the problem. However if I try to do something like


this gives out numbers like

{82.407, 131.036, -43.4576, -309.158, -510.697, -615.093, -660.662, -689.591, -725.895}

which are completely out of any possible result for that quantity (note that in this list appears again the quantity that I evaluated before, for $eFolds = 40$, and it is totally different). I chose to post the issue using Table instead than ParametricPlot, since the problems appear already using Table. Anyway, the code I tried to use also was


that gives, like table, very strange numbers, and

ParametricPlot[{nsplot[U3, c, n],rplot[U3, c, n]}, {c, 3, 5}, {n, 30, 70}]

that gives, too, a very strange plot.

To see if something was wrong with the code, I tried with a simpler potential (namely, $1-Cos[\frac{a}{b}]$) that can be solved "with pen and paper" and I tried to do the parametric plot firstly using the "analytical solution", and then using this code: the results are exactly identical.

Maybe the problem is in the small imaginary parts? If there's and imaginary part the $Cos$ in the potential will turn into some $exp$ and maybe that messes up the calculations? But then again if I simply calculate something like


it gives me the correct result. So I don't really know where's the problem.

Anyway, I hope you can give me some hint...



P.S.: I used "Compile->False" in NDSolve because I recalled that it was sometimes best to use it when dealing with numerical calculations, I hope that's not the problem...

I'll write down what I've found after experimenting a little more with the code, in the hope that someone wants to take a look.

The original code was a little different from the one I posted in the question (my fault, I really didn't imagine the problem would be in the part I left out... Well, lesson learned). The potential was written in this way

U[Fi_,qi_] = (-1)^{Fi+1}*Sum[Cos[n*qi*#1/#2]/n^5,{n,1,Infinity}];

so it was a "pure function" (I think they are called that way) only with respect to #1 and #2. I did that because in the future I would have needed to change those parameters Fi and qi a bit.

The rest of the notebook was exactly the same: the only differences were that I set those parameters Fi and qi to one for simplicity

phiendtemp[b_?NumericQ] := FindRoot[inthelp[U3,x]-Sqrt[2]*b,{x,48/10},MaxIterations->Infinity,Compiled->False][[1,2]]


phiendtemp[b_?NumericQ] := FindRoot[inthelp[U3[1,1],x]-Sqrt[2]*b,{x,48/10},MaxIterations->Infinity,Compiled->False][[1,2]]

and didn't give problems. The other quantity that had a different definition was


it too contained


instead that only


and it also didn't give any problems when I tried to use

 phicmbtemp2[b_,j_,eFolds_] = y[b,j][eFolds]/.phicmbtemp1;
 phicmb[b_,eFolds_] = phicmbtemp2[b,phiend[b],eFolds];

The rest was the same. Then when I tried to calculate the $x$ and $y$ coordinates I would write something like




the first command worked, the second gave the problems that I wrote in my first question. Then I tried to remove the dependence of the potential on those two parameters, thus making my notebook exactly like the one I posted in the question, and everything worked...

So my new questions are:

1) why was there a problem to begin with? I used other times (albeit with older versions of Mathematica) a construct like

U[Fi_,qi_] = (-1)^{Fi+1}*Sum[Cos[n*qi*#1/#2]/n^5,{n,1,Infinity}];

and never encountered problems;

2) how can I write a potential that uses some parameters in a way similar to what I wrote in point 1)? I'm starting to think that maybe I can define also those two parameters $Fi$ and $qi$ as #3 and #4... Maybe it is better that way. Anyways I'd really like to know why there are problems with that kind of construct.

share|improve this question
You have a few unbalanced parentheses in your first code block, and you need to group the first two equations with curly braces in your ParametricNDSolve. After fixing these errors, I can't reproduce your problem: evaluating nsplot[U3, 4, 40] gives me {1.00743 + 0. I} while evaluating the Table gives {{1.01355 + 0. I}, {1.00743 + 0. I}, {1.0038 + 0. I}, {1.00189 + 0. I}, {1.00093 + 0. I}}. Maybe you have some conflicting definitions left over; try restarting the kernel and evaluating the relevant bits again? – Rahul Apr 7 '14 at 23:54
Thanks for the answer. Regarding the curly braces in ParametricNDSolve, I forgot to add them in the question but on my notebook they are present. As far as the other braces are concerned my notebook doesn't tell me I miss braces so maybe I simply forgot to add them in the question as well (I'll check again anyways). I will try to restart the kernel but really I can't understand such a strange behaviour. – giova7_89 Apr 8 '14 at 15:16

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.