# Why Mathematica gives me different results depending on the order I set the parameters when I solve this differential equation?

I am asking Mathematica to solve a differential equation using DSolve, with two different parameters, q and m, and then evaluate at a point.

So I write the following lines of code:

x[y_] = w[y] /.
First[DSolve[{w'[y] == (1 - (2 m y^2)/(y^2 + q^2)^(3/2))^-1},
w[y], y]] /. C[1] -> 0;
m = 1.; q = 1/2;

x[2]


Where the output is (I just chose a random point to evaluate this function, it does not really matter where):

-4.60489 - 1.14725*10^-19i .


However, If I first set the parameters (using the same values as before) and then let Mathematica solve it:

m = 1.; q = 1/2;

x[y_] = w[y] /.
First[DSolve[{w'[y] == (1 - (2 m y^2)/(y^2 + q^2)^(3/2))^-1}, w[y],
y]] /. C[1] -> 0

x[2]


It yields:

Inactive[Integrate][(
0.25 Sqrt[1.  + 4. K[1]^2] +
1. K[1]^2 Sqrt[1.  + 4. K[1]^2])/(-4. K[1]^2 +
0.25 Sqrt[1.  + 4. K[1]^2] +
1. K[1]^2 Sqrt[1.  + 4. K[1]^2]), {K[1], 1, 2}]


which is an integral that cannot be analytically solved.

How come this happens? Why would the order that I define my constants change the outcome of how Mathematica handles this differential equation? What I find even weirder, is that (at least to me) it seems like it would be simpler to set the parameters first and then solve it.

• When inexact numbers appear in symbolic expression, Mathematica typically braches to numerical routines more suitable for numeric input. Apr 6, 2020 at 16:55
• Mathematica 101: Do not mix inexact numbers with solvers meant to give exact results (like DSolve, Integrate, Solve, etc...). Apr 6, 2020 at 22:02
• I think a more serious problem is, which result is correct? Or both of them are correct/incorrect? NIntegrate seems to suggest the last integral doesn't converge. Apr 7, 2020 at 3:19
• Note that the integrand is singular at four points on the real axis for the parameters chosen. So, the solutions obtained in the two cases may not be for the same regions of y. Also, C[1] is not necessarily the same quantity in the two solutions. Apr 7, 2020 at 13:53

Both solutions are correct, and the behavior of Mathematica is (to me) not entirely surprising. To begin, plot the right side of the ODE.

ComplexPlot3D[(1 - (2 m y^2)/(y^2 + q^2)^(3/2))^-1 /. {m -> 1, q -> 1/2}, {y, 3},
AxesLabel -> {Re[y], Im[y], None}, ImageSize -> Large, LabelStyle -> {15, Black, Bold},
PlotRange -> {0, 6}]


Four poles, on the real axis, and two branch cuts, on the imaginary axis, are visible. On the other hand, for m -> 1/2, the poles move off the real axis.

With m and q unspecified, as in the first case in the Question, Mathematica cannot know the location of the poles or where in the complex plane the desired values of y lie. So, it appears to make unknown assumptions and proceed. With m and q specified, it knows were the poles are but not on which side of the branch cuts the desired values of y lie. Again, it makes assumptions and proceeds. I do find it surprising, however, that Mathematica makes different assumptions for the second case in the Question and the case considered in Ulrich Neumann's Answer, which differ only by inserting m -> 1. and m -> 1.

Now, address which solutions are correct.

s1f = DSolveValue[w'[y] == (1 - (2 m y^2)/(y^2 + q^2)^(3/2))^-1, w, y];
Plot[Evaluate[{w'[y] /. w -> s1f, (1 - (2 m y^2)/(y^2 + q^2)^(3/2))^-1}
/. {m -> 1, q -> 1/2} // Chop], {y, -5, 5}, ImageSize -> Large,
PlotStyle->{Automatic, Dashed}, LabelStyle -> {15, Black, Bold}]


Visibly, the left and right sides of the ODE are equal. (Comparing the two with FullSimplify would be preferable, but I got tired of waiting for it to complete.) Doing the same for the second case in the Question and for Ulrich Neumann's Answer show that they too satisfy the ODE.

• You can add PlotStyle->{Automatic, Dashed} to Plot to make the identity easier to see :) . Apr 8, 2020 at 2:44
• @xzczd Good suggestion. Thanks. Apr 8, 2020 at 4:31

It seems to be dangerous to use y one both sides of the function definition x[y].

Perhaps this approach is what you're looking for

W[m_, q_] :=DSolveValue[{w'[y] == (1 - (2 m y^2)/(y^2 + q^2)^(3/2))^-1}, w , y]
W[1,1/2][2]
(*(-3.21859 + 0. I) + C[1]*)

• The function definition isn't relevant. If you use W[1., 1/2][2], the Inactive[Integrate][…] will be generated again. Apr 7, 2020 at 3:08