MathGroup Archive 2006

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

Search the Archive

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  
you get the same answers.

Andrzej Kozlowski


  • Prev by Date: Re: critical points of a third order polynomial fit (simplification)
  • Next by Date: Re: Re: Finding the Number of Pythagorean Triples below a bound
  • Previous by thread: Re: critical points of a third order polynomial fit (simplification)
  • Next by thread: Re: critical points of a third order polynomial fit (simplification)