Re: Root finding needs higher accuracy
- To: mathgroup at smc.vnet.net
- Subject: [mg123185] Re: Root finding needs higher accuracy
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 25 Nov 2011 04:58:40 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
----- Original Message -----
> From: "HwB" <hwborchers at googlemail.com>
> To: mathgroup at smc.vnet.net
> Sent: Thursday, November 24, 2011 5:57:06 AM
> Subject: Root finding needs higher accuracy
>
> 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
One possibility is to assess order of vanishing at the root. One can do this numerically as below. First we get an approximation. Machine arithmetic is actually okay for this, but I use higher precision just to make the zeros more apparent.
In[65]:= froot =
x /. FindRoot[f12[x] == 0, {x, 3/2}, WorkingPrecision -> 40,
AccuracyGoal -> Infinity, PrecisionGoal -> 20, MaxIterations -> 500]
Out[65]= 1.648721270699964636815140169116380771581
Define the derivative function.
In[66]:= fderiv[x_, n_] := D[f12[x], {x, n}]
Evaluate the first several. I include the zeroth (the function itself) so one will observe that the residual is zero to as many places as we used precision in the computations.
In[67]:= Table[fderiv[x, n], {n, 0, 5}] /. x -> froot
Out[67]= {0.*10^-40,
5.965503326769*10^-27, -7.29680399262098108582527420*10^-14, \
0.446260320296992429926744093547958318234, \
-0.812011699419998272665014156460529212842, \
1.970039966974547966294755465552043162439}
We see that this function very likely vanishes to third order at the root. That is, both first and second derivatives also vanish there. So we'll find that same root, for the second derivative. This is now well conditioned.
In[68]:= fderivroot =
x /. FindRoot[fderiv[x, 2] == 0, {x, 3/2}, WorkingPrecision -> 40,
AccuracyGoal -> Infinity, PrecisionGoal -> 30, MaxIterations -> 500]
Out[68]= 1.648721270700128146848650787814163571654
Since we happen to know the exact result, let's check against that.
In[69]:= fderivroot - Sqrt[E]
Out[69]= 0.*10^-40
Notice that we can also actually solve exactly for the first or second derivative vanishing.
In[71]:= Solve[fderiv[x, 1] == 0, x]
Out[71]= {{x -> Sqrt[E]}, {x -> Sqrt[E]}}
We could use this result to show that f12[x] also vanishes there, thus confirming exactly what we learned numerically.
Daniel Lichtblau
Wolfram Research