Re: Re: Re: FindRoot accuracy/precision
- To: mathgroup@smc.vnet.net
- Subject: [mg11510] Re: [mg11460] Re: Re: FindRoot accuracy/precision
- From: "C. Woll" <carlw@u.washington.edu>
- Date: Fri, 13 Mar 1998 12:22:00 -0500
Hi Mike,
The problem you are running across is that unless you are very lucky, it
is usually impossible to solve a set of 7 nonlinear equations
symbolically. Since Solve always attempts to solve problems
symbolically, it is unlikely that Mathematica will ever produce a
result for you. Instead, you should try using FindRoot as described by
Daniel Lichtblau in his response.
Carl Woll
Dept of Physics
U of Washington
On Thu, 12 Mar 1998, Michael Chang wrote:
> Thanks again for your response!
>
> Daniel Lichtblau (danl@wolfram.com) writes:
> > Michael Chang wrote:
> >> ...
> >> Hi! Thanks for your response!
> >>
> >> With regards to an example, suppose I have the following (made up on the
> >> spur of the moment).
> >>
> >> 100==x1*x2-x3*x4;
> >> 1*10^-3==x3/x2-x1^2;
> >> 2*10^-6==x4/x1+x2*x3;
> >> 5*10^-13==x1*x4-x3/x2;
> >>
> >> eq1={100,1*10^-3,2*10^-6,5*10^-13};
> >> eq2={x1*x2-x3*x4, etc.}
> >>
> >> (Note the great "spread" in numerical values.) I have four equations
> >> with four unknowns, and eq1 and eq2 *both* have *infinite* precision.
> >> Because I don't know where to start my parameter space search, I try
> >> using a random starting point. I'm looking for positive real roots
> >> since the 7 nonlinear equations that I'm *actually* dealing with come
> >> about from an electrical R & C circuit.
> >>
> >> Hence, I try the following:
> >>
> >> ok=0;While[ok==0, ans=FindRoot[(eq1-eq2)==0,{x1,Random[Real,x1max]},
> >> {x2,Random[Real,x2max]},{x3,Random[Real,x3max]},
> >> {x4,Random[Real,x4max]},MaxIterations->1000];
> >> If[Min[{x1,x2,x3,x4}/.ans]>=0,
> >> ok=1; Print[ans]]]
> >>
> >> Now since eq1 and eq2 both have infinite precision, should I be using
> >>
> >> Rationalize[Random[Real,ximax],0]
> >>
> >> in the above expressions? The main problem is that the actual value of
> >>
> >> (eq1-eq2)
> >>
> >> is "small" (the largest residual is of order 10^-2), but not small in
> >> percentage for the smaller valued elements. My default for
> >> WorkingPrecision is 16, but this doesn't seem to give (eq1-eq2) to be
> >> order 10^-16 ! I will try increasing WorkingPrecision, but am finding
> >> that the calculations are becoming increasingly slower.
> >>
> >> Once again, any suggestions and advice would be greatly appreciated!
> >>
> >> Thanks again,
> >>
> >> Mike
> >
> > Here are a few pointers that may help.
> >
> > (1) Often you can use Solve or NSolve. For example, try out exprs =
> > {x1*x2-x3*x4-100,
> > x3/x2-x1^2-10^(-3),
> > x4/x1+x2*x3-2*10^(-6),
> > x1*x4-x3/x2-5*10^(-13)};
> > vars = {x1,x2,x3,x4};
> > Timing[sol = Solve[exprs==0, vars];] nsol1 = N[sol, 100];
> > exprs /. nsol1 // Chop
> >
> > (2) The solutions to the above system are not real. Hence if we try
> > FindRoot and give real starting values, we will not converge to a
> > solution (this is because FindRoot uses some flavor of Newton
> > iterations, and all the data is real).
> >
> > If I do, say,
> >
> > max = 50;
> > start := Table[Random[Complex,max+I*max,100], {4}]
> > FindRoot[Evaluate[Thread[exprs==0]],
> > Evaluate[Sequence @@ Thread[{vars,start}]],
> > WorkingPrecision->100, AccuracyGoal->20,
> > MaxIterations->300]
> >
> > then I get a fairly decent result.
> >
> > (3) You do not need to rationalize the starting values. Indeed, the
> > precision of the starting values should be largely immaterial, assuming
> > the code is robust (I have not looked into this issue).
> >
> > (4) Given the spread of inputs, I would be surprised if machine
> > arithmetic will suffice for your needs. Hence setting WorkingPresision
> > to some bignum level is probably the appropriate thing to do. You will
> > most likely also want to up the values for AccuracyGoal and
> > MaxIterations, as I have done above. As you realize, this will result
> > in slower computation, but I see no easy remedy for that.
> >
> >
> > Daniel Lichtblau
> > Wolfram Research
>
> Okay, I've had Solve running for two days now trying to solve my actual
> seven variable problem, but haven't gotten any output yet. (I'm
> running a NEC P220MHz PC with Windoze 95 and 64MB RAM and Mathematica
> 3.0.)
>
> Specifically, my seven equations are as follows:
>
> eqs={369664/981 - 1/(x2*(x3+x4))==0,
> -6513786578938559/9223372036854775808+x1*x4+x1*x5+x6*x7==0,
> -3858337547701687/9444732965739290427392+x1*x2*x4*x5+x1*x4*x6*x7+
> x1*x5*x6*x7==0,
> -1/97280000000+x1*x2*x4*x5*x6*x7==0, 981/369664-x2*x3-x2*x4==0,
> 77/924160000-x1*x2*x3*x4-x1*x2*x3*x5-x2*x3*x4*x6-x2*x3*x6*x7-
> x2*x4*x6*x7==0,
> 1/1848320000000-x1*x2*x3*x4*x6*x7-x1*x2*x3*x5*x6*x7==0}
>
> (The fractions above may seem ridiculous, but they give me infinite
> precision.) So ... I try
>
> ans=Solve[eqs,{x1,x2,x3,x4,x5,x6,x7}]
>
> Thus far, I haven't had any output (unfortunately) after about 1.5 days.
> Now I'm a newbie at Mathematica, so is there something that I'm doing
> wrong here? I initially started just wanted to find positive real root
> solutions to these seven equations, but now I'm just struggling to find
> *any* complex solution to this set of equations. Sigh ...
>
> Again, I'd be grateful for any suggestions or help that anyone may have!
>
> Thanks again!
>
> Mike
>
>
>