8
$\begingroup$

I've read here that Mathematica 10 can Obtain symbolic solutions to delay differential equations.

Would that help in numericaly evaluating Dickman's function $\rho(u)$ ?

It is a delay differential equation defined for $u\in\mathbb R^+$ by $$\begin{cases} \rho(u)=1 & \text{if }0\le u\le1\\ u\cdot\rho'(u)+\rho(u-1)=0& \text{if }u>1 \end{cases}$$


Otherwise said, can Mathematica 10 (or something) help improve/simplify this (working) code?

Dickm$pr=40 (*desired ABSOLUTE (rather than relative) accuracy*);

Dickm$ex=Ceiling[Dickm$pr*Log[3,10]]+1 (*order of expansions; empirical*);

Dickm$ma=Ceiling[Dickm$ex/Log[Dickm$ex]] (*domain limit*);

(*Dickm$dv[n] approximates DickmanRho[t+n-1/2] for -1/2\[LessEqual]t\[LessEqual]1/2*)

Dickm$dv[1]=N[1,Dickm$pr+1];

Dickm$dv[n_]:=Dickm$dv[n]=Module[{u=Dickm$dv[n-1],v},v=Integrate[(u+O[\.b2]^Ceiling[Dickm$ex-n])/(\.b2+n-1/2),\.b2];(Normal[u]/.\.b2->1/2)+(Normal[v]/.\.b2->-1/2)-v];

DickmanRho[y_?NumberQ]:=Which[y<=0,1,y>Dickm$ma,N[1,y Log[10.,y]]-1,True,Normal[Dickm$dv[Ceiling[y]]]/.\.b2->y-Ceiling[y]+1/2];

(* Example use *)
DickmanRho[15]
$\endgroup$
1

1 Answer 1

7
$\begingroup$

Yes!

fun = NDSolve[{u ρ'[u] + ρ[u - 1] == 0, ρ[u /; u < 1] == 1}, ρ, {u, 0, 20}, 
    WorkingPrecision -> 40, PrecisionGoal -> ∞, AccuracyGoal -> 40, 
    Method -> "StiffnessSwitching"][[1, 1, 2]];

fun[15]
(* 7.5899080042980595046528227779709126741*10^-20 *)

DickmanRho[15]
(* 7.5899080042980595047*10^-20 *)

It works not only in V10 (in V8 at least).

$\endgroup$
3
  • 2
    $\begingroup$ @fgrieu, I tune the options for NDSolve and now it works much better. However, it takes several seconds which is ~100 times slower than your clever method! $\endgroup$
    – ybeltukov
    Sep 12, 2014 at 7:41
  • $\begingroup$ Looks like Method made a difference. $\;$ I wonder if we could make use the new DSolve capabilities. $\endgroup$
    – fgrieu
    Sep 12, 2014 at 8:06
  • 1
    $\begingroup$ @fgrieu I tried DSolve, but it produces very huge formulas so we have to manually simplify and rounding them. But it is what you have already done... $\endgroup$
    – ybeltukov
    Sep 12, 2014 at 8:14

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.