Services & ResourcesWolfram ForumsMathGroup Archive

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

ImplicitSolve

• To: mathgroup at smc.vnet.net
• Subject: [mg25760] ImplicitSolve
• From: Carl Woll <carlw at u.washington.edu>
• Date: Sat, 21 Oct 2000 14:43:13 -0400 (EDT)
• References: <200010190835.EAA21542@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

In this post, I want to describe a function which can solve a set of N
implicit equations for N dependent variables, where the implicit equations
also contain a single unknown parameter, which I will call the independent
variable. I have mentioned how one could do this in the past, but I hadn't
provided a function to actually carry out the procedure. This new function
is called ImplicitSolve, and has the usage message

?ImplicitSolve

ImplicitSolve[eqns, {x->x0,y->y0,...}, {x, xmin, xmax},
   opts] finds a solution to the implicit equations eqns
   for the functions y, ... with the independent
   variable x in the range xmin to xmax. The root
   {x0,y0,...} should satisfy the equations, or should
   provide a good starting point for finding a solution
   when using FindRoot. Currently, the only available
   option is AccuracyGoal, but a better ImplicitSolve
   would include the possibility of supplying options
   for both the FindRoot and NDSolve function calls.

I will give the function definition at the end of this post.

As an example, recently, Anesh Sooklal wanted to use Solve to
parametrically determine w as a function of k, where w depends implicitly
on k. This implicit dependence involved essentially a quartic equation,  so
Mathematica should be able to perform this task. See below

Anesh Sooklal wrote:

> I have been trying to solve the following equation using Solve,
> When I type in the following line I get a message saying that an illegal
> operation has been performed.
>
> f[w_]: = 1./(w^2)  + (.9/30.)/(w - 2. Sqrt[200./10.]k)^2   - 1.-
> (3./30.)/( Sqrt[200./10.]k^2 )
>
> Solve[f[w]==0., w]

It turns out that a rather serious bug occurs when the above code is
executed. The interesting thing here is that if all the machine numbers are
replaced with integers,  the Solve works fine, although the resulting
solution is rather messy, since the solution of a quartic is quite nasty.

Now, rather than using Solve, one could use the function ImplicitSolve.
First, since the equation has 4 roots for a typical value of k, we will
need to pick one of these roots as a starting point. Arbitrarily choosing a
root with k=1, we find

Solve[(f[w] /. k -> 1) == 0, w]

{{w -> -0.989151}, {w -> 0.989233}, {w -> 8.77187},

  {w -> 9.11659}}

Now that we have a root, we can use ImplicitSolve as follows:

ImplicitSolve[{f[w]==0},{k->1,w->.989233},{k,1,2}]

{{w -> InterpolatingFunction[{{1., 2.}}, <>]}}

What is returned is an interpolating function which gives the value of w
for any value of k between 1 and 2.

Carl Woll
Physics Dept
U of Washington

Here is the definition of ImplicitSolve:

Clear[ImplicitSolve]

Options[ImplicitSolve]={AccuracyGoal->6};

ImplicitSolve[eqn_List,rt_,{x_,min_,max_},opts___?OptionQ]:=Module[{x0,y,y0,accuracy},

    (* options *)
    accuracy=AccuracyGoal/.{opts}/.Options[ImplicitSolve];

    (* root *)
    x0=Cases[rt,(x->a_)->a];
    If[x0=={},
      Message[ImplicitSolve::badroot,x];Return[],
      x0=x0[[1]]];
    {y,y0}=Transpose[DeleteCases[rt,x->x0]/.Rule->List];

    (* check root *)
    If[!(And@@Thread[Abs[eqn/.Equal->Subtract/.rt]<10^-accuracy]),
     Message[ImplicitSolve::inaccurate];
     y0=FindRoot[
     Evaluate[eqn/.x->x0],
       Evaluate[Sequence@@(Transpose[{y,y0}])],opts][[All,2]],
      Null,
      Message[ImplicitSolve::incomplete];Return[]
    ];

    (* get interpolating function *)
    NDSolve[
     Flatten[{D[eqn/.Thread[Rule[y,#[x]&/@y]],x],Thread[#[x0]&/@y==y0]}],
     y,
     {x,min,max},
     PrecisionGoal->accuracy]
    ]

ImplicitSolve::usage="ImplicitSolve[eqns, {x->x0,y->y0,...}, {x, xmin,
xmax}, opts] finds a solution to the implicit equations eqns for the
functions y, ... with the independent variable x in the range xmin to xmax.
The root {x0,y0,...} should satisfy the equations, or should provide a good
starting point for finding a solution when using FindRoot. Currently, the
only available option is AccuracyGoal, but a better ImplicitSolve would
include the possibility of supplying options for both the FindRoot and
NDSolve function calls.";

ImplicitSolve::badroot="Supplied root is missing value for `1`";
ImplicitSolve::incomplete="Supplied root is incomplete";
ImplicitSolve::inaccurate="Supplied root is inaccurate, using FindRoot to
improve accuracy";

• References:
  • Solve...
    • From: Anesh Sooklal <anesh@fermi.udw.ac.za>