11
$\begingroup$

Bug introduced in 13.1 or earlier and persisting through 13.1.0 or later

I have a first-order ODE that I want to solve:

$$\frac{dy}{dt}=\frac{1}{2t}+\frac{y^2}{2}$$.

The result I get back from

DSolve[y'[t] == 1/(2 t) + 1/2 (y[t])^2, y, t]

is, in TeX form,

$$\frac{c_1 J_1(\sqrt{t})+\frac{\sqrt{t}}{2}\left[c_1 \left(J_0(\sqrt{t})-J_2(\sqrt{t})\right)-2J_0(\sqrt{t})\right]}{t \left(1-c_1\right)J_1(\sqrt{t})}$$

This includes the arbitrary constant of integration $c_1$, as is right and proper for a first-order ODE.

However, using FullSimplify, or indeed the recurrence relation for Bessel functions $x(J_{\nu+1}(x)+J_{\nu-1}(x))=2\nu J_{\nu}$, this can be reduced to

$$-\frac{J_0(\sqrt{t})}{\sqrt{t}\ J_1(\sqrt{t})}$$

Somehow the constant $c_1$ has disappeared from the problem. What this appears to be saying is that the solution to my ODE is somehow "rigid", in the sense that I can't solve it from an arbitrary initial condition. This contradicts my entire intuition about differential equations!

This leads me to think that what DSolve returns may not be the general solution to a differential equation.

So my question is: what does DSolve actually return? If it's not the general solution, then how might I find the general solution for my problem?

If it does return the general solution, then it would appear that my problem is pathological in some sense. In what sense though? Where can I find out more about this?

Improve this question
$\endgroup$
4
  • $\begingroup$ The plot thickens! Changing variables to $x=\sqrt{t}$, the ODE becomes $$\frac{dy}{dx}=\frac{1}{x}+x y^2$$ DSolve then gives the answer (after a FullSimplify): $$y=-\frac{1}{x}\frac{Y_0(x)+J_0(x)c_1}{Y_1(x)+J_1(x)c_1}$$ Not only does $c_1$ remain here (as I believe it should), but I also can't find a way to make this equivalent to the solution in terms of $\sqrt{t}$, which it should since they are linked by a simple change of variable! This system is strange. Any help would be very welcome! $\endgroup$
    – jms547
    Oct 13, 2022 at 20:32
  • $\begingroup$ I am now thinking the solution is correct, but FullSimplify is the one who made mistake because the solution actually verifies the ode according to Mathematica when substituted back into the ode. Any way, I give below derivation which gives solution that do not cancel the constant when simplified. If you report this to WRI they will know if it is FullSimplify the cause or not. $\endgroup$
    – Nasser
    Oct 13, 2022 at 22:28
  • $\begingroup$ The bug centers on the power $t^{-1}$: DSolve[{y'[t] == b/t + c (y[t])^2}, y[t], t] // FullSimplify has the same problem. -- Actually, the problem arises for negative powers of $t$. $\endgroup$
    – Michael E2
    Oct 14, 2022 at 13:36
  • $\begingroup$ Yet another workaround: y[t] -> ReplaceAll[t -> t - 1]@ DSolveValue[{y'[t] == 1/(2 (1 + t)) + 1/2 (y[t])^2}, y[t], t] // FullSimplify -- basically it seems these workarounds kick the execution path out of the buggy Ricatti solver. Please report the bug to WRI. $\endgroup$
    – Michael E2
    Oct 14, 2022 at 13:45

3 Answers 3

Reset to default
6
$\begingroup$

This is how to obtain the solution from the original $y(t)$ ode.

This is standard Riccati ode. Using known transformation (I can show the steps if needed) transforms the first order ode in $y(t)$ to second order ode in $u(t)$, which Mathematica solves correctly. Then the solution is transformed back to $y(t)$ which now gives a solution that keeps the constant of integration when fully simplified.

The transformation comes out to be
$$ y = -2 u'/u $$

This is the code

ClearAll[y, u, t];
odeInY = y'[t] == 1/(2 t) + 1/2 (y[t])^2; (*original ode *)
odeInU = odeInY /. y -> Function[{t}, -2*u'[t]/u[t]] // Simplify (*transformation*)

Mathematica graphics

solve the second order $u$ ode

solInU = u[t] /. First@DSolve[odeInU, u[t], t]

Mathematica graphics

Go back to $y(t)$

solInY = -2*D[solInU, t]/solInU // FullSimplify

Mathematica graphics

Divide numerator and denominator by $c_2$ and label $\frac{c_1}{c_2}$ as $c_0$ gives the solution in $y(t)$

solInY = solInY /. {C[2] -> 1, C[1] -> C[1]/C[2]} /. C[1]/C[2] -> C[0] //   FullSimplify

Mathematica graphics

Verify the above solution

odeInY /. y -> Function[{t}, Evaluate[solInY]] // FullSimplify

Mathematica graphics

Compare the above with what Mathematica gives

Which actually also verifies the ode also!

sol = DSolve[y'[t] == 1/(2 t) + 1/2 (y[t])^2, y, t]
odeInY /. sol // FullSimplify

Mathematica graphics

So it is possible that the solution is correct to the ode as written (just different form) but it is FullSimplify who is at fault!

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Check out DSolveChangeVariables[Inactive[DSolve][{y'[t] == 1/(2 t) + 1/2 (y[t])^2}, y, t], u, t, y[t] == -2 u'[t]/u[t]] // Simplify. (+1) $\endgroup$
    – Michael E2
    Oct 13, 2022 at 23:30
  • $\begingroup$ Thanks Nasser, this is a very helpful reply. I will just point out that I don't think FullSimplify is at fault as I can analytically reduce the first solution (in terms of J0 and J1 and with a C1) into its form without the C1. Using the recurrence relation for Bessel functions a factor of (1+C1) cancels from the top and bottom. So the problem is in DSolve somewhere. $\endgroup$
    – jms547
    Oct 15, 2022 at 16:32
  • $\begingroup$ Sure, I was not sure myself if it is FullSimplify or DSolve. WRI will know for sure when they look at it. Any way, the constant of integration should not go away ofcourse when simplified. $\endgroup$
    – Nasser
    Oct 15, 2022 at 18:30
5
$\begingroup$

The general solution should contain the homogeneous solution. Therefore consider the homogeneous equation:

y0 = y /. DSolve[y'[t] - 1/2 y[t]^2 == 0, y, t][[1]]

enter image description here

This is not contained in the result from MMA. Therefore the solution:

enter image description here

is not the general solution, but is a particular solution. This is also in line with the fact that c1 gets canceled in the result.

I think this is a bug and would report it to [email protected]

Improve this answer
$\endgroup$
3
  • $\begingroup$ OPs equation is not a linear differential equation. The terminology that you use (homogeneous equation, particular solution, …) is for linear differential equations and should not be used here. $\endgroup$
    – user293787
    Oct 14, 2022 at 2:41
  • $\begingroup$ @user293787 A homogeneous DE need not be linear. $\endgroup$
    – Daniel Huber
    Oct 14, 2022 at 10:38
  • $\begingroup$ This ODE may also be called separable. Homogeneous usually means invariant under $t\mapsto\lambda t,\,y\mapsto\lambda y$ positive or at least nonzero). This is one is invariant under $t\mapsto\lambda t,\,y\mapsto\lambda^{-1} y$, and such an ODE is often called quasihomogeneous with exponents $(1,-1)$. The change of variables $u=1/y$ converts it to a homogeneous ODE. $\endgroup$
    – Michael E2
    Oct 15, 2022 at 12:31
4
$\begingroup$

Daniel Huber's answer pointed out to me another path to a valid solution:

DSolveChangeVariables[
 Inactive[DSolve][{y'[t] == 1/(2 t) + 1/2 (y[t])^2}, y, t], u, t, 
 y[t] == 1/u[t]]
1/u[t] /. First@Activate@% // FullSimplify
(*
Inactive[DSolve][
 {-(Derivative[1][u][t]/u[t]^2) == 1/(2 t) + 1/(2 u[t]^2)}, 
 u, t]

(-2 BesselY[0, Sqrt[t]] - BesselJ[0, Sqrt[t]] C[1])/
  (Sqrt[t] (2 BesselY[1, Sqrt[t]] + BesselJ[1, Sqrt[t]] C[1]))
*)

I also agree with Daniel that this is a bug.

Updated to show the solution instead of it satisfying the ode

Another workaround is to disable the bug:

ode = {y'[t] == 1/(2 t) + 1/2 (y[t])^2};
Block[{DSolve`DSolveFirstOrderODEDump`Riccati = $Failed &},
 y[t] /. DSolve[ode, y, t] // FullSimplify
 ]
(*
{(-2 I BesselY[0, Sqrt[t]] - BesselJ[0, Sqrt[t]] C[1])/(
 Sqrt[t] (2 I BesselY[1, Sqrt[t]] + BesselJ[1, Sqrt[t]] C[1]))}
*)

This satisfies the ode, which can be checked by plugging it in and simplifying.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Maybe I misunderstand your edit, but after Quit[], also ode={y'[t]==1/(2 t)+1/2 y[t]^2}; ode/.DSolve[ode,y,t]//FullSimplify returns {{True}}, at least in V12.3. $\endgroup$
    – user293787
    Oct 14, 2022 at 14:02
  • $\begingroup$ @user293787 Apparently you didn't inspect the solution? $\endgroup$
    – Michael E2
    Oct 14, 2022 at 14:42
  • $\begingroup$ You are right, I did not, I was just reading your post. I think it is clearer after your latest update. $\endgroup$
    – user293787
    Oct 14, 2022 at 14:55
  • $\begingroup$ @user293787 I was originally going to post the y[x] formula, but since it differed from the other solution in form, I changed my mind and posted the test that the result was a solution. Your first comment made me realize one couldn't tell whether DSolve returned a nontrivial answer. THanks, $\endgroup$
    – Michael E2
    Oct 14, 2022 at 16:42
  • $\begingroup$ Thanks Michael. I'm using Mathematica 11 which doesn't have DSolveChangeVariables, but doing the variable change by hand I get the same result as you. Anyway your edit indeed suggests that the Riccati solver is the source of the bug, and I will report it. Thank you! $\endgroup$
    – jms547
    Oct 15, 2022 at 16:36

Your Answer

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

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