MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Root finding needs higher accuracy

  • To: mathgroup at smc.vnet.net
  • Subject: [mg123176] Re: Root finding needs higher accuracy
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Fri, 25 Nov 2011 04:57:02 -0500 (EST)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • References: <201111241157.GAA29024@smc.vnet.net>
  • Reply-to: drmajorbob at yahoo.com

Here's the function and an attempt using FindRoot:

Clear[f]
f[x_] = Log[x] + x^2/(2 Exp[1]) - 2 x/Sqrt[Exp[1]] + 1;
x /. Quiet@
    FindRoot[f[x], {x, 1.0, 3.4}, Method -> "Brent",
     WorkingPrecision -> 100, PrecisionGoal -> 25] // Timing
N[Sqrt@Exp@1 - Last@%, 25]

{0.003579, 1.64872127073570073780869588498038437600}

-3.557259096004509716622080*10^-11

Looking at the plot:

Plot[{f@x, f'[x], f''[x], f'''[x]}, {x, 1, 2}]

suggests that it's a double, even TRIPLE root, and if that's true, the  
solution is easy:

root = x /. First@Solve[f'[x] == 0, x]
f@root

Sqrt[E]

0

We can also solve this:

Solve[f''[x] == 0, x]

{{x -> -Sqrt[E]}, {x -> Sqrt[E]}}

Sqrt[E] is the root.

BUT, if we don't know it's a double root or Solve can't find a root of f',  
THEN what?

Here's a method that works for this function (and maybe a few others):

Clear[trisect]
trisect[{low_, high_}] := Module[{tbl, split},
   tbl = Table[{f@x,
      x}, {x, {low, low + (high - low)/3, low + 2 (high - low)/3,
       high}}];
   split = SplitBy[tbl, Negative];
   {split[[1, -1, -1]], split[[-1, 1, -1]]}
   ]

precision = 100;
iterations = 80;
Timing@Mean@
   Block[{$MinPrecision = precision, $MaxPrecision = precision},
    NestWhile[trisect, SetPrecision[{1, 2}, precision],
     Abs[Subtract @@ #] > 10^-25 &, 1]]

{0.012396, \
1.64872127070012814684865080152405879098770003155385410176215449001051\
1663472116782925272061059040521}

Choosing adequate "precision" is a problem in itself, of course.

Bobby

On Thu, 24 Nov 2011 05:57:06 -0600, HwB <hwborchers at googlemail.com> wrote:

> I would like to numerically find the root of the following function
> with up to 20 digits.
>
>     f12[x_] := Log[x] + x^2 / (2 Exp[1]) - 2 x / Sqrt[Exp[1]] + 1
>
> This problem is known to be difficult for solvers in double precision
> arithmetics. I thought it should be easy with Mathematica, but the
> following attempts were not successful.
>
>     SetPrecision[
>         x /. FindRoot[f12[x], {x, 1.0, 3.4}, Method -> "Brent",
>                  AccuracyGoal -> Infinity, PrecisionGoal -> 20], 16]
>     # 1.648732212532746
>     SetPrecision[
>         x /. FindRoot[f12[x], {x, 1.0, 3.4}, Method -> "Secant",
>                  AccuracyGoal -> Infinity, PrecisionGoal -> 20], 16]
>     # 1.648710202030051
>
> The true root obviously is Sqrt[Exp[1]]//N = 1.648721270700128...
>
> The symbolic solver explicitely says it cannot solve this expression.
> What do I need to do to get a much more exact result out of
> Mathematica?
>
> Many thanks, Hans Werner
>


-- 
DrMajorBob at yahoo.com



  • Prev by Date: Re: Using Equal with Real Numbers
  • Next by Date: Re: Non trivial substitution in a very long output
  • Previous by thread: Re: Root finding needs higher accuracy
  • Next by thread: Re: Root finding needs higher accuracy