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$
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).

| improve this answer | |
$\endgroup$
  • 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 '14 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 '14 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 '14 at 8:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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