Message Boards Message Boards

Unexpected result from DSolve? (Version dependence)

POSTED BY: Isao Maruyama
7 Replies

I see. However, I have never seen before this kind of incorrect numerical answer in Mathematica. For example, maxima's desolve has the auto-rationalize trick.

enter image description here

I thought Mathematica did it in a more elegant way than maxima. If DSolve requires exact numbers only, then I have to rewrite my teaching materials...

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
POSTED BY: Isao Maruyama

The original ODE is equivalent to a homogeneous 3rd-order ODE with a repeated root:

inflated = {D[2.3*r''[t] == 1 - r'[t], t],
   2.3*r''[t] == 1 - 0, r[0] == 0, r'[0] == 0}
(*  {2.3` r'''[t] == -r''[t], 2.3` r'[t] == 1, r[0] == 0, r'[0] == 0}  *)

First@inflated /. r -> Function[t, Exp[L t]] /. t -> 0 // Solve
(*  {{L -> -0.434783}, {L -> 0.}, {L -> 0.}}  *)

A repeated root is numerically unstable — any round-off error is likely to split the root in two. This is what happens in the procedure used by DSolve to solve the OP's nonhomogeneous 2nd-order ODE:

TracePrint[
 DSolve[{2.3*r''[t] == 1 - r'[t], r[0] == 0, r'[0] == 0} 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.80144*10^16 E^(5.55112*10^-17 t) - 2.3 E^(-0.434783 t) C[1] + C[2]}},
 {r[0] == 0, Derivative[1][r][0] == 0}, {C[1], C[2]}, {t}]

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

The number 5.55112*10^-17 in the first exponential term is a round-off error of 2.^-54 and 1.80144*10^16 is its reciprocal. If C[2] were -1.80144*10^16, they the two terms together would be approximately t (for small enough t). That would make this solution equivalent to the correct, exact solution.

The following shows the steps from the general to the particular solution, verifying the some of the description above:

gensol = DSolve[{23/10.*r''[t] == 1 - r'[t]}, r[t], t]
gensolFN = DSolve`DSolveToPureFunction[gensol]
(*
{{r[t] -> 
   1.80144*10^16 E^(5.55112*10^-17 t) - 2.3 E^(-0.434783 t) C[1] + C[2]}}

{{r -> Function[{t}, 
    1.80144*10^16 E^(5.55112*10^-17 t) - 2.3 E^(-0.434783 t) C[1] + C[2]]}}
*)

csol = Solve[{r[0] == 0, r'[0] == 0} /. First@gensolFN, {C[1], C[2]}]
gensol /. First@csol // Factor
% // Simplify
(*
{{C[1] -> -1., C[2] -> -1.80144*10^16}}

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

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

Round-off error and its effects is normally considered one of the risks of working with floating-point numbers. As others have pointed out, DSolve[] is an exact solver, and it has (probably) not been designed to be careful with approximate floats. OTOH, it solves the higher order system without trouble:

DSolve[inflated, r[t], t] // Simplify
(*  {{r[t] -> -2.3 + 2.3 E^(-0.434783 t) + 1. t}}  *)

Evidently, it's using a different sequence of steps, but one should expect that, given one is homogeneous and the other nonhomogeneous. You could report it to Wolfram support. They might be able to improve DSolve[].

POSTED BY: Michael Rogers

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 6 days 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
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract