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