0
$\begingroup$

I've got the following piece of code:

    ClearAll[Evaluate[Context[] <> "*"]];

    f[x_] := N[-E^x ExpIntegralEi[-N[x]]];
    g[x_] := NIntegrate[f[y]/y^2, {y, N[x], ∞}];
    addzero[r_, a_] := f[1/a] - f[1/r];
    lb[r_, a_, γ_] := (g[1/a] - g[1/r])/(r + γ);

    diff[r_, γ_] := N[lb[r, γ + r, γ] - addzero[r, γ + r]];
    ropt[γ_] := FindRoot[diff[r, γ], {r, γ/4, γ/2}][[1]][[2]];

    Off[NIntegrate::inumr];
    Off[NIntegrate::ncvb];
    Off[NIntegrate::slwcon];

    For[γ = 10, γ <= 500, γ += 10^
        Floor[Log10[γ]],
         Print["γ=", γ, ", \!\(\*SubscriptBox[\(r\), \(opt\)]\)=", ropt[γ];];
    ];

The code basically solves a certain nonlinear equation. The equation involves a special function (viz. the exponential integral) plus an integral of that special function over an unbounded region (so there is a risk that it might be numerically unstable). For the curious, it's not difficult to see from the code what the equation is. The code works but Mathematica begins to complain when one of the parameters ($\gamma$) of the equation becomes sufficiently high so that the integral of the special function is problematic to evaluate numerically.

I am not a Mathematica guru, but I suspect Mathematica can potentially offer higher accuracy. Any suggestions as to how to tweak the code so as to make it capable of handling larger values of $\gamma$ (say, on the order of $10^3\div 10^4$)?

I am using Mathematica 10.0.1.0, 64bit, Student edition.

Thanks in advance.

$\endgroup$
2
$\begingroup$

The easiest solution is to give the root finder [FindRoot] better starting points.

For example use:

ropt[γ_] := FindRoot[diff[r, γ], {r, 2, 4}][[1]][[2]];

Also your print statement has a minor syntax error.

For[γ = 10, γ <= 5000, γ += 10^Floor[Log10[γ]], 
Print["γ=", γ, ", \!\(\*SubscriptBox[\(r\), \(opt\)]\)=", ropt[γ]];];

(* γ=50000, Subscript[r, opt]=2.29546 *)
$\endgroup$
  • $\begingroup$ Thanks much, that seems to work. And thanks for pointing out the syntax error, I corrected it in the question. $\endgroup$ – Alex Dec 6 '14 at 19:01

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.