MathGroup Archive 2008

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

Search the Archive

Re: Finding a continuous solution of a cubic

  • To: mathgroup at smc.vnet.net
  • Subject: [mg85941] Re: [mg85899] Finding a continuous solution of a cubic
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Thu, 28 Feb 2008 02:51:07 -0500 (EST)
  • References: <200802270925.EAA16050@smc.vnet.net> <0B1C17D7-E5DC-44DE-BB5F-94F043A34A66@mimuw.edu.pl>

On 27 Feb 2008, at 18:36, Andrzej Kozlowski wrote:

>
> On 27 Feb 2008, at 10:25, Hugh Goyder wrote:
>
>> I am investigating the maxima and minima of the expression, e, below.
>> The expression is a quartic in x and has two parameters a and b  
>> which,
>> for my problem, are restricted to the range -1 < a < 1 and 0 < b (the
>> interesting part is 0<b<1).
>>
>> I start with a simple approach and take the derivative with respect  
>> to
>> x and solve to find the three turning points. The solution is
>> complicated and not easy to examine. I work out the discriminant to
>> see when I will have three real roots or one and plot the region. I
>> work out three functions for the three roots r1, r2 and r3, and three
>> functions for the value of e at the roots, f1, f2, and f3. I am
>> expecting three real functions when within the region defined by the
>> discriminant and one outside this region. I plot the three roots and
>> the three functions (there are some artifacts where the functions  
>> turn
>> from real to complex, this is not a problem).
>>
>> Now my problem. The function f3 is continuous in the region where
>> there are three solutions and this I like. Functions f1 and f2 are  
>> not
>> continuous but contain holes. If functions f1 and f2 are plotted
>> together then they each fill in the others holes. The holes makes my
>> life difficult because I do not have continuous functions.
>>
>> Is it possible to make Solve find roots that will give rise to three
>> functions one corresponding to the maxima and two corresponding to
>> each minima without the functions crossing over and jumping between
>> the turning points?
>>
>> I also notice, from the plots when shown together,  that my desired
>> functions are even with respect to the parameter a. Is this correct
>> and how would I show this if I can't assemble them as continuous
>> functions?
>>
>> Thanks
>>
>> Hugh Goyder
>>
>> e = ((1 - x)^2 + b*(-1 + a)^2)*((1 + x)^2 +
>>         b*(1 + a)^2);
>>
>> d = Simplify[D[e, x]];
>>
>> dis = Discriminant[d, x];
>>
>> ContourPlot[dis, {a, -1, 1}, {b, 0, 1},
>>  ContourShading -> False, Contours -> {0},
>>  FrameLabel -> {"a", "b"}]
>>
>> sol = Solve[d == 0, x];
>>
>> ClearAll[x1, x2, x3, f1, f2, f3];
>> x1[a_, b_] := Evaluate[x /. sol[[1]]];
>> x2[a_, b_] := Evaluate[x /. sol[[2]]];
>> x3[a_, b_] := Evaluate[x /. sol[[3]]];
>> f1[a_, b_] := Evaluate[e /. sol[[1]]];
>> f2[a_, b_] := Evaluate[e /. sol[[2]]];
>> f3[a_, b_] := Evaluate[e /. sol[[3]]];
>>
>> r1 = Plot3D[x1[a, b], {a, -1, 1}, {b, 0, 1}]
>>
>> r2 = Plot3D[x2[a, b], {a, -1, 1}, {b, 0, 1}]
>>
>> r3 = Plot3D[x3[a, b], {a, -1, 1}, {b, 0, 1}]
>>
>> Show[r1, r2, r3, PlotRange -> All]
>>
>> s1 = Plot3D[f1[a, b], {a, -1, 1}, {b, 0, 1}]
>>
>> s2 = Plot3D[f2[a, b], {a, -1, 1}, {b, 0, 1}]
>>
>> s3 = Plot3D[f3[a, b], {a, -1, 1}, {b, 0, 2}]
>>
>> Show[s1, s2]
>>
>>
>>
>
> Mathematica will certainly not give you such a formula. I think in  
> your case it will be quite hard to obtian one, and I can't at the  
> moment suggest a good way to do so, but it may be helpful to first  
> consider the the cubic geometrically and understand how Mathematica  
> deals with the ordering of its roots. In general a cubic (with real  
> coefficients) can have either 1 or all 3 real roots.  If the cubic  
> ir represented by a polynomial poly (in your case poly = d) in the  
> variable x  then these roots can are given by the Mathematica  
> functions Root[poly,x,1], Root[poly,x,2] and Root[poly,x,3].
> For example:
>
> In[31]:= Root[x^3 + x + 1, x, 1]
> Out[31]= Root[#1^3 + #1 + 1 &, 1]
>
> In[32]:= ToRadicals[%]
> Out[32]= -(2/(3*(-9 + Sqrt[93])))^(1/3) + ((1/2)*(-9 +  
> Sqrt[93]))^(1/3)/3^(2/3)
>
>
>
> Root[poly,x,1],Root[poly,x,2] and Root[poly,x,3] are your functions  
> f1,f2, and f3.
>
> What Mathematica counts as the first root of a cubic, as mentioned  
> above, is always  real. So suppose you define your real root  
> function as Root[poly,x,1] and also suppose that we "start" in a  
> situation (a set of parameter values) where we have  one real root  
> and two conjugate non-real ones. The real root is our Root[poly,x, 
> 1]. Now, however, as we vary the parameters, the two non-real roots  
> Root[poly,x,2] and Root[poly,x,3] will at some point coalesce to  
> form a double real root, which may be smaller than your original  
> root. At this point the roots are "renamed", in other words these  
> two roots become Root[poly,x,1] and Root[poly,x,2] and the original  
> root then becomes Root[poly,x,3]. So if we stick to the function  
> Root[poly,x,1] we will have a discontinuous jump. In order to avoid  
> it, we should take Root[poly,x,3] as the value of our function after  
> the jump point.
>
> So in fact, to construct a continuous function you only need to  
> consider Root[poly,x,1] and Rot[poly,x,3] for the root will never  
> jump from the first to second place, because whenever two complex  
> roots coalesce, two new real of the same magnitude roots are created  
> and so what used to be root number one will become root number 3.
>
> As you noted, all the jumps take place on the curve  
> dis=Discriminant[d, x] = 0, because the choice of a smallest root  
> can only be changed once a double root appears. However, not every  
> appearance of a double root will result in a jump, because the two  
> real roots that appear due to the merging of two complex ones may be  
> larger than the real root we had before their merger.
>
> Now, in the case of a cubic with a single parameter, the equation  
> dis can actually be solved numerically so we can find all the jump  
> points and use the function Piecewise to combine them to get a  
> continuous function of the parameter. This is easy and has been  
> considered more than once on this list in the past. Your two  
> parameter case seems much more complicated (however, I have not  
> considered it very carefully!) since the discriminant is now a real  
> curve of degree 9 in the parameter space so at this time I have no  
> good idea to suggest. Still, I hope you can find something useful in  
> the above remarks.
>
> Andrzej Kozlowski


Just one more thing. You can see graphically how the "jumping" takes  
place by plotting together just the first and the third root of the  
cubic. Here is one way to do it:

F[a_, b_] :=
  Module[{poly, x}, poly = 4 x^3 + 4 a^2 b x + 4 b x - 4 x - 8 a b;
   Piecewise[{{{Root[poly, x, 1], Root[poly, x, 1]},
      Im[Root[poly, x, 3]] != 0}}, {Root[poly, x, 1], Root[poly, x,  
3]}]]

Plot3D[F[a, b], {a, 0, 1}, {b, -1, 1}]

Of course the vertical part should not be there. The "single" sheeted  
surface is the part given by Root[poly,x,1], which I had to tell  
Mathematica do draw twice, in order to have two values for every pair  
{a,b}. The double sheeted part of the surface consists of Root[poly,x, 
1] and Root[poly,x,3] over the part of the parameter space where there  
are three real roots. The upper part, that joins continusly onto the  
single sheeted part is Root[poly,x,3]. The lower part which ought to  
be isolated from the rest of the graph (but is to it connected by  
vertical lines) is Root[poly,x,1], so you can see where one has to  
change the choice of root to obtain a continuous graph.

Andrzej Kozlowski 


  • Prev by Date: NDSolve of a 3rd Order Nonlinear equation
  • Next by Date: Re: Series solution for complicated function
  • Previous by thread: Re: Finding a continuous solution of a cubic
  • Next by thread: v6: Input Cell Word Wrapping: how to set global default?