Re: critical points of a third order polynomial fit (simplification)

• To: mathgroup at smc.vnet.net
• Subject: [mg68489] Re: [mg68466] critical points of a third order polynomial fit (simplification)
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Mon, 7 Aug 2006 01:41:10 -0400 (EDT)
• References: <200608060656.CAA23460@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```On 6 Aug 2006, at 08:56, Chris Chiasson wrote:

> Here are some commands that find the (two) critical locations of a
> cubic polynomial fitted to four points. The four points are given as
> bracket={y1,x1,y2,x2,y3,x3,y4,x4}, where y is the range and x is the
> domain location. The code produces two sets of critical points. Both
> solution sets involve substituting the coefficients obtained from the
> four point "fit" into the general solution of the critical points
> (from calculus D[poly,x]==0). The only difference is that one of them
> calls a Simplify command that the other doesn't (compare rep[3] with
> rep[4]).
>
> Why are the answers for the second critical point so different? Do you
> think the simplification made in rep[3] is invalid? Does your computer
> give the same results that mine gives?
>
> In[1]:=
> \$Version
>
> Out[1]=
> 5.2 for Microsoft Windows (June 20, 2005)
>
> In[2]:=
> bracket={0.040911501042171394,9.797734083340343,0.003643048928782312,
>       9.939642325021731,0.0007478278717134088,10.027346441664564,
>       0.028647147834538127,10.169254683345951};
>
> In[3]:=
> eqn[1]=Plus@@Table[Times[c[i],Power[x,i-1]],{i,4}]\[Equal]y;
> eqn[2]=eqn[1]/.{x\[Rule]x[#],y\[Rule]y[#]}&/@Range[4];
> rep[1]=FullSimplify@Solve[eqn[2],Array[c,{4}]];
> rep[2]=FullSimplify@Solve[D[eqn[1],x],x];
> rep[3]=Simplify[rep[2]/.rep[1][[1]]];
> rep[4]=rep[2]/.rep[1][[1]];
> rep[5]={x[num_]\[RuleDelayed]bracket[[2 num]],
>       y[num_]\[RuleDelayed]bracket[[2 num-1]]};
>
> In[10]:=
> x/.rep[3]/.rep[5]//InputForm
>
> Out[10]//InputForm=
> {-1.4838127878952383*^10,
>  17.999998818597952}
>
> In[11]:=
> x/.rep[4]/.rep[5]//InputForm
>
> Out[11]//InputForm=
> {-1.4838127878952383*^10,
>  9.999999069105009}
>
> --
> http://chris.chiasson.name/
>

It's all caused the evil demon that often wakes up when you mix
symbolic algebra with machine precision arithmetic.

You have to be vary careful when using machine precision numbers with
any sort of symbolic algebra. In general, you cannot expect
algebraically equivalent expression to return the same answer when
machine precision numbers are substituted for variables. So just
replace your list of machine reals "bracket" by

bracketRat = Rationalize[bracket, 0];

define rep[6] by

rep[6] = {x[num_] :> bracketRat[[2 num]],
y[num_] :> bracketRat[[2 num - 1]]}

and

evaluate N[x /. rep[3] /. rep[6] // InputForm, 20]
N[x /. rep[3] /. rep[6] // InputForm, 20]

(you can use precision lower than 20 of course) and you will see that