Re: Re: FindRoot accuracy/precision
- To: mathgroup@smc.vnet.net
- Subject: [mg11364] Re: [mg11347] Re: FindRoot accuracy/precision
- From: Daniel Lichtblau <danl@wolfram.com>
- Date: Sat, 7 Mar 1998 02:06:37 -0500
- References: <199803030411.XAA02345@smc.vnet.net.> <6dk5kd$oam$3@dragonfly.wolfram.com> <199803060541.AAA22147@smc.vnet.net.>
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
- References:
- FindRoot accuracy/precision
- From: bt585@FreeNet.Carleton.CA (Michael Chang)
- Re: FindRoot accuracy/precision
- From: bt585@FreeNet.Carleton.CA (Michael Chang)
- FindRoot accuracy/precision