Group Abstract Group Abstract

Message Boards Message Boards

Unexpected result from DSolve? (Version dependence)

Could you help me to understand the following?

In[1]:= $VersionNumber

Out[1]= 14.

In[2]:= DSolve[{2.3*r''[t]==1-r'[t],r[0]==0,r'[0]==0},r[t],t]//FullSimplify

Out[2]= {{r[t] -> -2. + 2.3 E^(-0.434783 t)}}

In[3]:= DSolve[{m*r''[t]==f-r'[t],r[0]==0,r'[0]==0},r[t],t]/.{m->2.3,f->1}//FullSimplify

Out[3]= {{r[t] -> -2.3 + 2.3 E^(-0.434783 t) + t}}

2.3 seems to be a magic number. Is this a bug?

(2024-10-23 ADD) Thank you for replies. I found version dependence. I report the result in the following by using Gianluca Gorni's DiscretePlot.

derivativeAtZero[a_] := DSolveValue[{a*r''[t] == 2 - r'[t], r[0] == 0, r'[0] == 0}, r, t]'[0];
fncAtZero[a_] := DSolveValue[{a*r''[t] == 2 - r'[t], r[0] == 0, r'[0] == 0}, r, t][0];
DiscretePlot[{fncAtZero[a], derivativeAtZero[a]}, {a, Subdivide[1., 3, 200]}, PlotLabel -> $VersionNumber]

enter image description here

Correct answer(right figure) : version=13.1, 13.0, 12.2?

Wrong answer(left figure) : version=14.1, 14

(2024-10-23 ADD2) Thanks to Michael Rogers' reply, I noticed a way to reveal version-dependent numerical instability.

The equation 2.3 r’’(t)=1 - r’(t) can be generalized to 2.3 r’’(t)=exp(k t) - r’(t). The following command

DSolveValue[{2.3*r''[t]==Exp[k*t]-r'[t],r[0]==0,r'[0]==0},r[t],t]//FullSimplify//FullForm//Print

gives version-dependent outputs

version 14: Times[-1,Plus[0.4347826086956521`,Times[-0.43478260869565216`,Power[E,Times[k,t]]],Times[0.9999999999999999`,k],Times[-0.9999999999999999`,Power[E,Times[-0.43478260869565216`,t]],k]],Power[Plus[Times[0.43478260869565216`,k],Times[1.`,Power[k,2]]],-1]]
version 13: Times[-1,Plus[0.4347826086956522`,Times[-0.4347826086956522`,Power[E,Times[1.`,k,t]]],Times[1.`,k],Times[-1.`,Power[E,Times[-0.4347826086956522`,t]],k]],Power[Plus[Times[0.4347826086956522`,k],Power[k,2]],-1]]

One can see small deviations of the same number 0.4347..., that is 1/2.3. In version 14, the numerical instability 0.4347826086956521 - 0.43478260869565216 = 5.55112*10^-17 appears in the solution shown by Michael Rogers.

This numerical instability depends on the number 2.3. For example, we obtain the correct answer and no numerical instability for the number 2.31 even if we use version 14.

DSolveValue[{2.31*r''[t]==1-r'[t],r[0]==0,r'[0]==0},r[t],t]//FullSimplify
(*
-2.31 + 2.31 E^(-0.4329 t) + 1. t
*)
DSolveValue[{2.31*r''[t]==Exp[k*t]-r'[t],r[0]==0,r'[0]==0},r[t],t]//FullSimplify//FullForm//Print
(*
Times[-1,Plus[0.4329004329004329`,Times[-0.4329004329004329`,Power[E,Times[k,t]]],Times[1.`,k],Times[-1.`,Power[E,Times[-0.4329004329004329`,t]],k]],Power[Plus[Times[0.4329004329004329`,k],Power[k,2]],-1]]
*)

(2024-10-24 ADD) In conclusion, I have understood that this unexpected result from DSolve was due to the numerical instability of a round-off error existing only in latest versions(>=14) of Mathematica. In the following, I report a new magic number with large round-off error (worst case), aWorst, by showing version-dependent solutions of $a\times v'(t)=e^{k t} - v(t) $ with using the following commands.

aWorst=5750000000002614783/2500000000000000000;
eq[a_,k_]:=a*v'[t]==Exp[10^-k*t]-v[t];
check[a_]:=Block[{},
Print[TraditionalForm[eq[a,k]],"\n a=",FullForm[a]," with Version",$VersionNumber,"\n",
TableForm[{#,DSolveValue[eq[a,#],v[t],t]}&/@Range[20],TableHeadings->{None,{"k","v(t)"}}]]];
check[2.3]
check[N[aWorst]]
check[2.31]

The result for version 14 shows unexpected round-off error while that for version 13 is reasonable.

enter image description here enter image description here To avoid it, one can use the following quick hack.

dSolveValue[eqs_,fnc_,arg_]:=Expand[DSolveValue[Rationalize[eqs,0],fnc,arg]]/.n_?NumberQ:>N[n];

I guess something has happened during the version update, and it might be not intended by the developer. I will report this to Wolfram support. I would like to thank you for your help and replies.
POSTED BY: Isao Maruyama
7 Replies
Sort By:
Bill Nelson
Bill Nelson
Posted 1 year ago

Try changing 2.3 to 23/10 and remove the FullSimplify and see if you can correctly manually simplify the result and get the same thing that FullSimplify gives you.

Then try substituting the solution into each of your three original equations and use D to calculate each needed derivative and see if the result is True in every case.

Both those appear to work for me.
POSTED BY: Bill Nelson

Very strange, the result is totally discontinuous in the parameter:

derivativeAtZero[a_] := 
  DSolveValue[{a*r''[t] == 2 - r'[t], r[0] == 0, r'[0] == 0},
     r, t]'[0];
DiscretePlot[derivativeAtZero[a], {a, Subdivide[1., 3, 200]}]
POSTED BY: Gianluca Gorni

DSolve, Solve are symbolic solvers,works correct with only exact numbers.

Exact numbers have infinity precision.

    DSolve[{23/10*r''[t] == 1 - r'[t], r[0] == 0, r'[0] == 0}, r[t], 
      t] // Simplify
    DSolve[{m*r''[t] == f - r'[t], r[0] == 0, r'[0] == 0}, r[t], 
       t] /. {m -> 23/10, f -> 1} // Simplify

   (*{{r[t] -> -(23/10) + 23/10 E^(-10 t/23) + t}}*)
   (*{{r[t] -> -(23/10) + 23/10 E^(-10 t/23) + t}}*)

Regards M.I.
POSTED BY: Mariusz Iwaniuk
POSTED BY: Michael Rogers
POSTED BY: Isao Maruyama

Thanks. It works also for me. A quick hack is the following;

dSolve[eqs_, fnc_, arg_] := DSolve[Rationalize[eqs], fnc, arg];
dSolve[{2.3*r''[t] == 1 - r'[t], r[0] == 0, r'[0] == 0}, r[t], t]
POSTED BY: Isao Maruyama

Thank you for your detailed analysis about numerical instability.

As you mentioned, I could get the correct answer by using the higher order system in version 14.0, which is magic for me:

DSolve[{D[2.3*r''[t] == 1 - r'[t], t], 2.3*r''[0] == 1 - r'[0], r[0] == 0, r'[0] == 0}, r[t], t] // Simplify

And, I noticed that Mathematica ver.13.0 gives different (correct) answer:

TracePrint[ DSolve[{2.3*r''[t] == 1 - r'[t], r[0] == 0, r'[0] == 0}, r[t],   t], _DSolve`DSolveBoundaryValueProblem, TraceForward -> True,  TraceInternal -> True]
(*
DSolve`DSolveBoundaryValueProblem[DSolve,DSolve`DSolveParserDump`result,DSolve`DSolveParserDump`initial,DSolve`DSolveParserDump`ctab,DSolve`DSolveParserDump`independents]

DSolve`DSolveBoundaryValueProblem[DSolve,{{r[t]->1. t-2.3 E^(-0.434783 t) Subscript[\[ConstantC], 1]+Subscript[\[ConstantC], 2]}},{r[0]==0,(r^\[Prime])[0]==0},{Subscript[\[ConstantC], 1],Subscript[\[ConstantC], 2]},{t}]

{r[t] -> 
   1. E^(-0.434783 t) (2.3 - 2.3 E^(0.434783 t) + 
      1. E^(0.434783 t) t)}}
*)
POSTED BY: Isao Maruyama
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard