Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 1998

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

Search the Archive

Re: Re: FindRoot accuracy/precision

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
  • Prev by Date: Re: Re: Simplifying algebraic expr: howto?
  • Next by Date: Re: Simplifying algebraic expr: howto?
  • Prev by thread: Re: Re: FindRoot accuracy/precision
  • Next by thread: Re: Re: Re: FindRoot accuracy/precision