Re: FindRoot in 5.0
- To: mathgroup at smc.vnet.net
- Subject: [mg42861] Re: [mg42819] FindRoot in 5.0
- From: Eric Mockensturm <emm10 at psu.edu>
- Date: Wed, 30 Jul 2003 19:31:30 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Richard,
Thanks for the reply. In the real application (the code below was just
a simple example to illustrate the problem) I can get FindRoot to
return a result using ?NumericQ as you suggest. However, it takes
about twice as long in 5.0 to find the root as it did in 4.2. I did
not use the SetSystemOptions command.
The actual code is below (please don't chuckle at the crudeness of it
;-). In 5.0 it takes roughly 2.5 sec to find the root and in 4.2 it
takes roughly 0.8 sec. Any suggestions? Could 5.0 be slower and
solving the linear equations in the tempfunc?
Eric
--
Eric Mockensturm, Ph.D.
Assistant Professor
Department of Mechanical & Nuclear Engineering
The Pennsylvania State University
SetDelayed[func1[Pattern[\[Zeta]darg, Blank[]], Pattern[\[Kappa]arg,
Blank[]], Pattern[\[CapitalOmega]arg, Blank[]],
Pattern[\[Tau]aarg, Blank[]], Pattern[\[Mu]rarg, Blank[]],
Pattern[\[Zeta]arg, Blank[]], Pattern[\[Mu]sarg, Blank[]],
Pattern[\[Rho]arg, Blank[]]], CompoundExpression[Set[\[Zeta]d,
\[Zeta]darg], Set[\[Kappa]b, \[Kappa]arg],
Set[\[CapitalOmega], \[CapitalOmega]arg], Set[\[Tau]a, \[Tau]aarg],
Set[\[Mu]r, \[Mu]rarg], Set[\[Zeta], \[Zeta]arg],
Set[\[Mu]s, \[Mu]sarg], Set[\[Rho], \[Rho]arg], Clear[\[Phi], t1],
Set[A, List[List[0, 1, 0, 0, 0, 0], List[Times[-1, Plus[\[Kappa]b,
2]], 0, \[Kappa]b, 0, 0, 0], List[0, 0, 0, 1, 0, 0],
List[Times[\[Kappa]b, Power[\[Mu]r, -1]], 0, Times[-1,
Times[\[Kappa]b, Power[\[Mu]r, -1]]], 0, 0, 0], List[0, 0, 0, 0, 0, 1],
List[0, 0, 0, 0, 0, Times[-1, \[Zeta]]]]], Set[List[eval, evec],
N[Eigensystem[A]]],
Set[cd, List[cd1, cd2, cd3, cd4, cd5, cd6]], Set[CC[Pattern[t,
Blank[]]],
Table[Times[Part[cd, ii], Power[E, Times[Part[eval, ii], t]]],
List[ii, 1, 6]]],
If[Equal[\[Zeta], 0], Set[CC[Pattern[t, Blank[]]], Join[Part[CC[t],
Range[1, 5]], List[Times[cd6, t]]]]],
Set[\[Theta]phd[Pattern[t, Blank[]]], Sum[Times[Part[cd, ii],
Part[evec, ii, 1], Power[E, Times[Part[eval, ii], t]]],
List[ii, 1, 4]]], Set[\[Theta]rhd[Pattern[t, Blank[]]],
Sum[Times[Part[cd, ii], Part[evec, ii, 3], Power[E,
Times[Part[eval, ii], t]]], List[ii, 1, 4]]],
Set[\[Theta]shd[Pattern[t, Blank[]]], Plus[cd6, Times[cd5, t]]],
Set[\[Theta]ppd[Pattern[t, Blank[]]], Times[\[CapitalTheta]p,
Power[E, Times[I, Plus[\[Phi], Times[t, \[CapitalOmega]]]]]]],
Set[\[Theta]rpd[Pattern[t, Blank[]]], Times[\[CapitalTheta]r,
Power[E, Times[I, Plus[\[Phi], Times[t, \[CapitalOmega]]]]]]],
Set[\[Theta]spd[Pattern[t, Blank[]]], Times[\[CapitalTheta]s,
Power[t, 2]]],
Set[\[CapitalTheta]dsol,
Part[Solve[List[Simplify[Equal[Plus[Times[Plus[\[Kappa]b, 2],
\[CapitalTheta]p],
Times[-1, Times[\[Kappa]b, \[CapitalTheta]r]],
Times[\[CapitalTheta]p, Times[-1, Power[\[CapitalOmega], 2]]]],
\[Rho]]],
Simplify[Equal[Plus[Times[\[Kappa]b, Plus[\[CapitalTheta]r,
Times[-1, \[CapitalTheta]p]]],
Times[-1, Times[Power[\[CapitalOmega], 2], \[Mu]r,
\[CapitalTheta]r]]], 0]],
Equal[Times[2, \[Mu]s, \[CapitalTheta]s], Times[-1, \[Tau]a]]],
List[\[CapitalTheta]p, \[CapitalTheta]r, \[CapitalTheta]s]],
1]], Set[\[Theta]pd[Pattern[t, Blank[]]],
Chop[ReplaceAll[Plus[\[Theta]phd[t], \[Theta]ppd[t]],
\[CapitalTheta]dsol]]],
Set[\[Theta]rd[Pattern[t, Blank[]]],
Chop[ReplaceAll[Plus[\[Theta]rhd[t], \[Theta]rpd[t]],
\[CapitalTheta]dsol]]],
Set[\[Theta]sd[Pattern[t, Blank[]]],
Chop[ReplaceAll[Plus[\[Theta]shd[t], \[Theta]spd[t]],
\[CapitalTheta]dsol]]],
Set[Ae, List[List[0, 1, 0, 0], List[Times[-1, Plus[\[Kappa]b, 2]],
0, \[Kappa]b, 0], List[0, 0, 0, 1],
List[Times[\[Kappa]b, Power[Plus[\[Mu]r, \[Mu]s], -1]], 0,
Times[-1, Times[\[Kappa]b, Power[Plus[\[Mu]r, \[Mu]s], -1]]],
Times[-1, \[Zeta]]]]], Set[List[eval, evec], N[Eigensystem[Ae]]],
Set[ce, List[ce1, ce2, ce3, ce4, ce5, ce6]],
Set[\[Theta]phe[Pattern[t, Blank[]]], Chop[Sum[Times[Part[ce, ii],
Part[evec, ii, 1], Power[E, Times[Part[eval, ii], t]]],
List[ii, 1, 4]]]], Set[\[Theta]rhe[Pattern[t, Blank[]]],
Chop[Sum[Times[Part[ce, ii], Part[evec, ii, 3], Power[E,
Times[Part[eval, ii], t]]], List[ii, 1, 4]]]],
Set[\[Theta]ppe[Pattern[t, Blank[]]], Plus[Cp, Times[Power[E,
Times[I, Plus[\[Phi], Times[t, \[CapitalOmega]]]]],
\[CapitalTheta]p]]], Set[\[Theta]rpe[Pattern[t, Blank[]]],
Plus[Cr, Times[Power[E, Times[I, Plus[\[Phi], Times[t,
\[CapitalOmega]]]]], \[CapitalTheta]r]]],
Set[\[CapitalTheta]esol, Part[Solve[List[Equal[Expand[Plus[Times[I,
\[Zeta]d,
Plus[\[CapitalTheta]p, Times[-1, \[CapitalTheta]r]],
\[CapitalOmega]], Times[Plus[\[Kappa]b, 2], \[CapitalTheta]p],
Times[-1, Times[\[Kappa]b, \[CapitalTheta]r]],
Times[\[CapitalTheta]p, Times[-1, Power[\[CapitalOmega], 2]]]]],
\[Rho]],
Equal[Expand[Plus[Times[I, \[Zeta]d, Plus[\[CapitalTheta]r,
Times[-1, \[CapitalTheta]p]], \[CapitalOmega]],
Times[Plus[\[CapitalTheta]r, Times[-1, \[CapitalTheta]p]],
\[Kappa]b],
Times[\[CapitalTheta]r, Plus[\[Mu]r, \[Mu]s], Times[-1,
Power[\[CapitalOmega], 2]]]]], 0]],
List[\[CapitalTheta]r, \[CapitalTheta]p]], 1]],
Set[Cesol, Part[Solve[List[Equal[Expand[Plus[Times[Plus[\[Kappa]b,
2], Cp], Times[-1, Times[\[Kappa]b, Cr]]]], 0],
Equal[Expand[Times[\[Kappa]b, Plus[Cr, Times[-1, Cp]]]],
Times[-1, \[Tau]a]]], List[Cr, Cp]], 1]],
Set[\[Theta]pe[Pattern[t, Blank[]]],
ReplaceAll[ReplaceAll[Plus[\[Theta]phe[t], \[Theta]ppe[t]],
\[CapitalTheta]esol], Cesol]],
Set[\[Theta]re[Pattern[t, Blank[]]],
ReplaceAll[ReplaceAll[Plus[\[Theta]rhe[t], \[Theta]rpe[t]],
\[CapitalTheta]esol], Cesol]],
Set[\[Theta]se[Pattern[t, Blank[]]],
ReplaceAll[ReplaceAll[Plus[\[Theta]rhe[t], \[Theta]rpe[t]],
\[CapitalTheta]esol], Cesol]],
Set[cdsol,
Join[Chop[Simplify[Chop[Part[Solve[List[Equal[\[Theta]pe[t1],
\[Theta]pd[t1]],
Equal[Derivative[1][\[Theta]pe][t1],
Derivative[1][\[Theta]pd][t1]], Equal[\[Theta]re[t1], \[Theta]rd[t1]],
Equal[Derivative[1][\[Theta]re][t1],
Derivative[1][\[Theta]rd][t1]]], List[cd1, cd2, cd3, cd4]], 1]]]],
Simplify[Part[Solve[List[Equal[\[Theta]re[t1], \[Theta]sd[t1]],
Equal[Derivative[1][\[Theta]re][t1],
Derivative[1][\[Theta]sd][t1]]], List[cd5, cd6]], 1]]]],
Set[T, Times[2, Pi, Power[\[CapitalOmega], -1]]],
SetDelayed[tempfunc[PatternTest[Pattern[argt1, Blank[]], NumericQ],
PatternTest[Pattern[arg\[Phi], Blank[]], NumericQ]],
CompoundExpression[Set[t1, argt1], Set[\[Phi], arg\[Phi]],
Set[cesol,
Part[Solve[Chop[Simplify[N[ReplaceAll[List[Equal[\[Theta]pe[0],
\[Theta]pd[T]],
Equal[\[Theta]re[0], \[Theta]rd[T]],
Equal[Derivative[1][\[Theta]pe][0], Derivative[1][\[Theta]pd][T]],
Equal[Derivative[1][\[Theta]re][0],
Derivative[1][\[Theta]rd][T]]], cdsol]]]], List[ce1, ce2, ce3, ce4]],
1]],
List[Plus[Re[ReplaceAll[Derivative[1][\[Theta]re][0], cesol]],
Times[-1, Re[ReplaceAll[ReplaceAll[Derivative[1][\[Theta]sd][T],
cdsol], cesol]]]],
Re[ReplaceAll[Plus[\[Tau]a, Times[\[Zeta],
Derivative[1][\[Theta]re][t1]], Times[\[Mu]s,
Derivative[2][\[Theta]re][t1]]],
cesol]]]]], Null]]
func1[0, 0.0564*30, 0.20, 0.0564, 0.5, 0, 20, 1]
Timing[FindRoot[tempfunc[xx, yy] == 0, {xx, 16, 16.5}, {yy, 2, 3}]]
On Wednesday, July 30, 2003, at 09:38 AM, Richard Gass wrote:
>
> On Wednesday, July 30, 2003, at 04:07 AM, Eric Mockensturm wrote:
>
>> I realize that FindRoot has changed but there seems to be a
>> fundamental difference in it between pre-5.0 versions and 5.0 that
>> breaks most of my notebooks. Maybe I've been following bad
>> 'programming' practices all these years, but....
>>
>> Anyway, a simple illustration follows.
>>
>> In 4.2:
>>
>> In[7]:=func[x_]:=(y=x;NIntegrate[BesselY[1/2,s],{s,1/2,y}])
>> In[8]:=FindRoot[(Print[x];func[x]),{x,2,4}]
>> Out[8]={x\[Rule]3.04729}
>>
>> In 5.0:
>>
>> In[1]:=func[x_]:=(y=x;NIntegrate[BesselY[1/2,s],{s,1/2,y}])
>> In[2]:=FindRoot[(Print[x];func[x]),{x,2,4}]
>> Out[2]={x\[Rule]4.}
>>
>> with errors...
>>
>> NIntegrate::nlim: s = x is not a valid limit of integration
>>
>> It seems that 5.0 is trying to evaluate func[x] symbolically before
>> start to search for a root. How do I make a global change to this
>> behavior. Does it have something to do with the HoldAll attribute?
>> FindRoot was been HoldAll in 4.2, too.
>>
>>
>>
>>
>
> This effects more than just FindRoot but you can restore the old
> behavior by using
> Developer`SetSystemOptions["EvaluateNumericalFunctionArgument" ->
> False] . This restores the old behavior but you take a substantial
> speed hit. In your case you can fix the problem (without the speed
> hit) by using
>
>> func[x_?NumericQ]:=(y=x;NIntegrate[BesselY[1/2,s],{s,1/2,y}])
>
> and then
>> FindRoot[func[x],{x,2,4},EvaluationMonitor:>Print[x]]
> or better yet
> points={}
> FindRoot[func[x],{x,2,4},EvaluationMonitor:>AppendTo[points,x]]
> ListPlot[points]
>
>
> Richard Gass
> Department of Physics
> University of Cincinnati
>
>