Re: Re: 2 stage solution with FindRoot

*To*: mathgroup at smc.vnet.net*Subject*: [mg57628] Re: [mg57449] Re: 2 stage solution with FindRoot*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Thu, 2 Jun 2005 05:17:15 -0400 (EDT)*References*: <d76och$7pd$1@smc.vnet.net> <200505280939.FAA21706@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Mukhtar Bekkali wrote: > Correction by example: > > I have functions f1[x1,x2,y1,y2] and f2[x1,x2,y1,y2]. I need to solve > the problem the following way. First, I need to solve for > {x1,x2}={{x1,x2}|D[f1,x1]=0&&D[f2,x2]=0}. A solution is > {x1[y1,y2],x2[y1,y2]}. Substituting into original functions yields > f1[x1[y1,y2],x2[y1,y2],y1,y2] and f2[x1[y1,y2],x2[y1,y2],y1,y2]. Then, > I need to solve for > {y1,y2}={{y1,y2}|D[f1[x1[y1,y2],x2[y1,y2],y1,y2],y1]=0&&D[f2[x1[y1,y2],x2[y1,y2],y1,y2],y2]=0} > > I apoligize for misspecification of the problem. I'll take an example that is solvable algebraically in order to show how to approach this in one step via FindRoot. The latter method will apply even when you cannot solve explicitly for x's in terms of y's. We'll work with f1[x1_,x2_,y1_,y2_] := 3*x1^2*x2 + 5*x2*y1*y2 - x1*y2 + y1 - 4*x1 + 2; f2[x1_,x2_,y1_,y2_] := 2*x1*x2^2*y2 - x2*y1*y2 + 3*x2*y1 + 2*x1 - 3*x2^2 + 5*y2 - 9; For generality I'll set partials to constants other than zero, but the same idea will apply to your case. In[4]:= InputForm[s1 = Solve[ {D[f1[x1,x2,y1,y2],x1], D[f2[x1,x2,y1,y2],x2]} == {3,-2}, {x1,x2}]] Out[4]//InputForm= {{x1 -> (3*(7 + y2))/(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2), x2 -> (6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2)/18}} Here are the functions rewritten explicitly in terms of {y1,y2}. In[5]:= InputForm[newf = {f1[x1,x2,y1,y2], f2[x1,x2,y1,y2]} /. s1[[1]]] Out[5]//InputForm= {2 + y1 - (12*(7 + y2))/(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2) - (3*y2*(7 + y2))/(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2) + (3*(7 + y2)^2)/(2*(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2)) + (5*y1*y2*(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2))/18, -9 + 5*y2 + (6*(7 + y2))/(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2) + (y1*(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2))/6 - (y1*y2*(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2))/18 + (y2*(7 + y2)*(6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2))/54 - (6 + 9*y1 + 14*y2 - 3*y1*y2 + 2*y2^2)^2/108} I'll solve for {y1,y2} and work with the second solution (which I already checked is real valued). soln = Solve[ {D[newf[[1]],y1],D[newf[[2]],y2]} == {-1,4}, {y1,y2}] [[2]]; We'll put together the two parts to evaluate for {x1,x2,y1,y2}, then show the numeric approximation. In[7]:= N[fullsoln = {x1,x2,y1,y2} /. First[s1] /. soln] Out[7]= {-1.42337, -0.49 The full solution, in unsimplified form, is a bit large and contains several dozen Root[] objects, that is, algebraic numbers. In[8]:= fullsoln//LeafCount Out[8]= 11971 In[10]:= Count[fullsoln, Root, Infinity, Heads->True] Out[10]= 98 Okay, so how do we do this when we cannot do the first step of solving for {x1,x2} in terms of {y1,y2}? This is a multivariable calculus exercise in chasing down partial derivatives. Clear[f1,f2] We set up initial equations. These are expressions we regard as being set to zero. partialsoln = {D[f1[x1,x2,y1,y2],x1]-c1, D[f2[x1,x2,y1,y2],x2]-c2}; We now recast this as having {x1,x2} explicit functions of {y1,y2}. We need this in order to solve for paritals of {x1,x2} with respect to {y1,y2}. implicitsoln = partialsoln /. {x1->x1[y1,y2], x2->x2[y1,y2]}; We use these to find, symbolically, those partials of {x1,x2}. These will be expressed in terms of (partials of {f1,f2} with respect to the various arguments. y1derivatives = First[Solve[D[implicitsoln,y1]==0, {D[x1[y1,y2],y1],D[x2[y1,y2],y1]}]]; y2derivatives = First[Solve[D[implicitsoln,y2]==0, {D[x1[y1,y2],y2],D[x2[y1,y2],y2]}]]; yderivatives = Join[y1derivatives,y2derivatives]; Now we form the full set of four equations. In[32]:= InputForm[exprs = Join[partialsoln, {D[f1[x1[y1,y2],x2[y1,y2],y1,y2],y1]-d1, D[f2[x1[y1,y2],x2[y1,y2],y1,y2],y2]-d2} /. yderivatives /. {x1[y1,y2]->x1, x2[y1,y2]->x2}]] Out[32]//InputForm= {-c1 + Derivative[1, 0, 0, 0][f1][x1, x2, y1, y2], -c2 + Derivative[0, 1, 0, 0][f2][x1, x2, y1, y2], -d1 + Derivative[0, 0, 1, 0][f1][x1, x2, y1, y2] - (Derivative[0, 1, 0, 0][f1][x1, x2, y1, y2]* (Derivative[1, 0, 1, 0][f1][x1, x2, y1, y2]*Derivative[1, 1, 0, 0][f2][ x1, x2, y1, y2] - Derivative[0, 1, 1, 0][f2][x1, x2, y1, y2]* Derivative[2, 0, 0, 0][f1][x1, x2, y1, y2]))/ (Derivative[1, 1, 0, 0][f1][x1, x2, y1, y2]*Derivative[1, 1, 0, 0][f2][x1, x2, y1, y2] - Derivative[0, 2, 0, 0][f2][x1, x2, y1, y2]* Derivative[2, 0, 0, 0][f1][x1, x2, y1, y2]) - (Derivative[1, 0, 0, 0][f1][x1, x2, y1, y2]* (Derivative[0, 2, 0, 0][f2][x1, x2, y1, y2]*Derivative[1, 0, 1, 0][f1][ x1, x2, y1, y2] - Derivative[0, 1, 1, 0][f2][x1, x2, y1, y2]* Derivative[1, 1, 0, 0][f1][x1, x2, y1, y2]))/ (-(Derivative[1, 1, 0, 0][f1][x1, x2, y1, y2]*Derivative[1, 1, 0, 0][f2][ x1, x2, y1, y2]) + Derivative[0, 2, 0, 0][f2][x1, x2, y1, y2]* Derivative[2, 0, 0, 0][f1][x1, x2, y1, y2]), -d2 + Derivative[0, 0, 0, 1][f2][x1, x2, y1, y2] - (Derivative[0, 1, 0, 0][f2][x1, x2, y1, y2]* (Derivative[1, 0, 0, 1][f1][x1, x2, y1, y2]*Derivative[1, 1, 0, 0][f2][ x1, x2, y1, y2] - Derivative[0, 1, 0, 1][f2][x1, x2, y1, y2]* Derivative[2, 0, 0, 0][f1][x1, x2, y1, y2]))/ (Derivative[1, 1, 0, 0][f1][x1, x2, y1, y2]*Derivative[1, 1, 0, 0][f2][x1, x2, y1, y2] - Derivative[0, 2, 0, 0][f2][x1, x2, y1, y2]* Derivative[2, 0, 0, 0][f1][x1, x2, y1, y2]) - (Derivative[1, 0, 0, 0][f2][x1, x2, y1, y2]* (Derivative[0, 2, 0, 0][f2][x1, x2, y1, y2]*Derivative[1, 0, 0, 1][f1][ x1, x2, y1, y2] - Derivative[0, 1, 0, 1][f2][x1, x2, y1, y2]* Derivative[1, 1, 0, 0][f1][x1, x2, y1, y2]))/ (-(Derivative[1, 1, 0, 0][f1][x1, x2, y1, y2]*Derivative[1, 1, 0, 0][f2][ x1, x2, y1, y2]) + Derivative[0, 2, 0, 0][f2][x1, x2, y1, y2]* Derivative[2, 0, 0, 0][f1][x1, x2, y1, y2])} We are ready to solve numerically. We redefine the actual expressions for {f1,f2}, and plug in the values {c1,c2,d1,d2} to which we set the partial derivatives. f1[x1_,x2_,y1_,y2_] := 3*x1^2*x2 + 5*x2*y1*y2 - x1*y2 + y1 - 4*x1 + 2; f2[x1_,x2_,y1_,y2_] := 2*x1*x2^2*y2 - x2*y1*y2 + 3*x2*y1 + 2*x1 - 3*x2^2 + 5*y2 - 9; exprsnum = exprs /. {c1->3,c2->-2,d1->-1,d2->4}; We'll pick starting values reasonably close to the actual solution in order to discourage FindRoot from discovering one different from that we already found. In[39]:= FindRoot[exprsnum==0, {x1,-3/2}, {x2,-1/2}, {y1,1/2}, {y2,-2}] Out[39]= {x1 -> -1.42337, x2 -> -0.497428, y1 -> 0.488355, y2 -> -2.75184} So we successfully recovered the solutions we found algebraically in two solve steps. What is the upshot of all this? A numeric approach wll begin at partialsolns=... above. No strenuous algebraic solving is required; we only solved a linear system to find partials of {x1,x2} in terms of {y1,y2}. We formulated one system of four equations in four variables, corresponding to the given specifications (except I set the various partials to constants other than zero, but this works fine if you use zero instead). One FindRoot sufficed to solve for all four variables. As with any FindRoot-based solving, I would guess that how well this works will depend on quality of your initial values. Daniel Lichtblau Wolfram Research