Re: Differentiating Functions and Root objects [ was Re: ArcCos[]]
- To: mathgroup at smc.vnet.net
- Subject: [mg24887] Re: [mg24833] Differentiating Functions and Root objects [ was Re: ArcCos[]]
- From: "Allan Hayes" <hay at haystack.demon.co.uk>
- Date: Sat, 19 Aug 2000 04:46:00 -0400 (EDT)
- References: <shNm5.8859$Ub4.81356@ralph.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Andrzej Kozlowski" <andrzej at tuins.ac.jp> wrote in message news:shNm5.8859$Ub4.81356 at ralph.vnet.net... > Allan: > > 1. I am afraid that you allowed yourself to be tricked by FindRoot :) > [rest of message copied below] Thanks Andrzej -yes I fell for that one - should have been suspicious of the large value - still, a salutary warning. For the record, here are two simpler examples of this situation: exp = 1/(1 + x^4); FindRoot[exp, {x, 0.1}] {x -> 39.2833} exp2 = x^4 + .00000001; FindRoot[exp2, {x, 0.1}] {x -> 0.023615} This does raise the question of what to count as a root in the presence of uncertainty, and what to do about roots at infinity. -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay at haystack.demon.co.uk Voice: +44 (0)116 271 4198 Fax: +44 (0)870 164 0565 "Andrzej Kozlowski" <andrzej at tuins.ac.jp> wrote in message news:shNm5.8859$Ub4.81356 at ralph.vnet.net... > Allan: > > 1. I am afraid that you allowed yourself to be tricked by FindRoot :) > > The root you found {t -> -818.986} is a pure illusion. The only roots are > the ones given in my earlier message on the same topic: -1.01505 and > 1.01505. You can see this simply by plotting the function: > > In[44]:= > Plot[Evaluate[D[Root[-t + 2*#1 + 2*t^2*#1 + #1^3 &, 1],t]],{t,-10,10}] > > But actually Mathematica can prove it algebraicalluy and very fast too using > Experimental`CylindricalAlgebraicDecomposition. > > In[5]:= > f[t_]:=D[Root[-t + 2*#1 + 2*t^2*#1 + #1^3 &, 1],t] > In[6]:= > Experimental`CylindricalAlgebraicDecomposition[f[t]==0,t] > Out[6]= > 2 4 > t == Root[-1 - 32 #1 + 32 #1 & , 1] || > > 2 4 > t == Root[-1 - 32 #1 + 32 #1 & , 2] > In[7]:= > N[%] > Out[7]= > t == -1.01505 || t == 1.01505 > > > 2. In general, a Root object depending on parameters has a complicated > branching structure with respect to these parameters and at certain points > will not be continuous, and hence certainly not differentiable (as a complex > function). Other than that (and barring bugs) differentiation of root > objects wiht respect to parameters is valid. However note one pitfall. The > case considered here is rather special: we are looking at a cubic f[t_]:= > Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & with real coefficients, and > Root[f[t],1]--its real root. The function can be viewed either as a real > function or a complex function. In the first case case, the derivative > D[f[t],t] exists at 0 and we get: > > In[8]:= > D[f[t],t]/.t->0 > Out[8]= > 1 > - > 2 > > > However both > > In[9]:= > Limit[(f[t]-f[0])/t,t->0] > > and > > In[10]:= > f[t]+O[t]^2 > > fail, since the complex derivative of f at 0 does not exist. > -- > Andrzej Kozlowski > Toyama International University, JAPAN > > For Mathematica related links and resources try: > <http://www.sstreams.com/Mathematica/> > > > > on 8/15/00 9:04 AM, Allan Hayes at hay at haystack.demon.co.uk wrote: > > > John's posting stirred me into taking another look a this problem. > > Here are four observations > > > > 1. We can differentiating a pure function with respect to a parameter > > > > D[Function[x, Sin[y] x^2], y] > > > > Function[x, x^2*Cos[y]] > > > > > > D[Sin[y] #^2 &, y] > > > > Cos[y]*#1^2 & > > > > 2. The original problem, differentiating Root objects > > > > a = Root[-t + 2*#1 + 2*t^2*#1 + #1^3 &, 1] > > > > Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1] > > > > (b = D[a, t]) // InputForm > > > > (1 - 4*t*Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1])/ > > (2 + 2*t^2 + 3*Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1]^ > > 2) > > > > Make this output explicit > > > > (rb = ToRadicals[b]) // InputForm > > > > (1 - 4*t*(-(2^(1/3)*(6 + 6*t^2))/ > > (3*(27*t + Sqrt[729*t^2 + 4*(6 + 6*t^2)^3])^ > > (1/3)) + (27*t + Sqrt[729*t^2 + 4*(6 + 6*t^2)^3])^ > > (1/3)/(3*2^(1/3))))/(2 + 2*t^2 + > > 3*(-(2^(1/3)*(6 + 6*t^2))/ > > (3*(27*t + Sqrt[729*t^2 + 4*(6 + 6*t^2)^3])^ > > (1/3)) + > > (27*t + Sqrt[729*t^2 + 4*(6 + 6*t^2)^3])^(1/3)/ > > (3*2^(1/3)))^2) > > > > Maybe Solve and NSolve will work on rb == 0, -- I gave up after a a brief > > wait -- but we can quickly get > > > > FindRoot[rb == 0, {t, 0}] > > > > {t -> -818.986} > > > > Is b correct for the rate of change of a respect to t? > > We can check > > > > rb - D[ToRadicals[a], t] // Simplify > > > > 0 > > > > We might expect with sol the corresponding polynomial would have a multiple > > root. > > > > Check: > > > > (poly = -t + 2*#1 + 2*t^2*#1 + #1^3 &[x] /. sol[[1]]) // InputForm > > > > 818.985893312185 + 1.3414777868887153*^6*x + x^3 > > > > > > Solve[% == 0, x] // InputForm > > > > {{x -> -0.0006105102159100438}, > > {x -> 0.0003052551079550219 - 1158.2218211072502*I}, > > {x -> 0.0003052551079550219 + 1158.2218211072502*I}} > > > > > > 3. Can we still get the solution without using ToRadicals? > > > > sol2 = FindRoot[b == 0, {t, 0}] > > > > FindRoot::"jsing": "Encountered a singular Jacobian at the point \!\(t\) > > = \!\ > > \(0.`\). Try perturbing the initial point(s)." > > > > FindRoot::"jsing": "Encountered a singular Jacobian at the point \!\(t\) > > = \!\ > > \(0.`\). Try perturbing the initial point(s)." > > > > FindRoot[b == 0, {t, 0}] > > > > However > > > > sol2 = FindRoot[b == 0, {t, 0.1}] > > > > {t -> -1.0150511440294088} > > > > A different solution. > > > > And we get > > > > Solve[(-t + 2*#1 + 2*t^2*#1 + #1^3 &[x] /. sol2[[1]]) == 0, x] > > > > {{x -> -0.2462928577523092}, > > {x -> 0.1231464288761546 - 2.0263644239932934*I}, > > {x -> 0.1231464288761546 + 2.0263644239932934*I}} > > > > 4. Is Differentiation of Root objects always valid? Is this method a general > > one for finding which parameter values give multiple roots? > > > > Allan > > --------------------- > > Allan Hayes > > Mathematica Training and Consulting > > Leicester UK > > www.haystack.demon.co.uk > > hay at haystack.demon.co.uk > > Voice: +44 (0)116 271 4198 > > Fax: +44 (0)870 164 0565 > > > > "John D. Hendrickson" <jdh at hend.net> wrote in message > > news:8n7q9r$1j7 at smc.vnet.net... > >> I've got Mathematica 4.0.2.0. Since I'm running Win95 I don't mind > >> crashing - thats pretty much a constant, so I tried your dare:) But mine > >> didn't crash - it gave me output. Also - I have a substitute that might > > be > >> what you want. > >> > >> In[1]:= > >> a = Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1]; > >> > >> In[2]:= > >> b = D[a, t]; > >> > >> In[4]:= > >> InputForm[ Solve[b == 0, t] ] > >> > >> Solve::"tdep": "The equations appear to involve the variables to be solved > > \ > >> for in an essentially non-algebraic way." > >> > >> Out[4]//InputForm= > >> Solve[(1 - 4*t*Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1])/ > >> (2 + 2*t^2 + 3*Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1]^ > >> 2) == 0, t] > >> > >> =================================================================== > >> Would the following be acceptible in your circumstance? > >> I assumed by making the unexplicit root object explicit: > >> ==================================================================== > >> > >> In[3]:= > >> Exit[] > >> > >> In[1]:= > >> a = Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1]; > >> > >> In[2]:= > >> ax = First[a][x] > >> > >> Out[2]= > >> \!\(\(-t\) + 2\ x + 2\ t\^2\ x + x\^3\) > >> > >> In[3]:= > >> b = D[ax, t]; > >> > >> In[4]:= > >> InputForm[ Solve[b == 0, t] ] > >> > >> Out[4]//InputForm= > >> {{t -> 1/(4*x)}} > >> > >> =================================================================== > >> > >> Allan Hayes wrote in message <8mtc0m$7qb at smc.vnet.net>... > >>> > >>> > >>> > >>> "Gianluca Gorni" <gorni at dimi.uniud.it> wrote in message > >>> news:8mqvd7$18j at smc.vnet.net... > >>> > >>> > >>> ----------- Snip ------------- > >>>> An unrelated problem: the following instructions consistently crash > >>>> my Mac Mathematica 4 kernel: > >>>> > >>>> a = Root[-t + 2*#1 + 2*t^2*#1 + #1^3 & , 1]; > >>>> b = D[a, t]; > >>>> Solve[b == 0, t] > >>>> > >>> -- > >>> Gianluca , > >>> It crashes Mathematica 4.02 on Windows also. > >>> > >>> Allan > >>> --------------------- > >>> Allan Hayes > >>> Mathematica Training and Consulting > >>> Leicester UK > >>> www.haystack.demon.co.uk > >>> hay at haystack.demon.co.uk > >>> Voice: +44 (0)116 271 4198 > >>> Fax: +44 (0)870 164 0565 > >>> > >>> > >>> > >>> > >> > >> > >> > > > > > > > > > > > > > > > > > > >