ComplexRoots
- To: mathgroup at yoda.physics.unc.edu
- Subject: ComplexRoots
- From: akyildiz at maths.ox.ac.uk ( Yilmaz Akyildiz tel 2-73558)
- Date: Thu, 19 Nov 92 21:04:08 GMT
(*
This is a Mathematics V.2.0 programme allowing the user to
numerically compute all the common roots of the two real functions
f1 and f2 of two variables x1 and x2 on a given rectangle in the
(x1, x2)-plane. It can easily be generalized to higher dimensions.
comroots[{f1, f2}, {x1, a, b}, {x2, c, d}, tolerance]
will find all the common roots of the functions f1(x1, x2) and
f2(x1, x2) with accuracy correct to the tolerance. It uses a version
of Interval Bisection Method.
For example:
comroots[{x1 Sin[x2] -1, x2 Cos[x1] + 1},
{x1, 0, 5}, {x2, 0, 9}]
finds all the roots of {x1 Sin[x2] -1 = 0, x2 Cos[x1] + 1 = 0}
in the intervals 0 <= x1 <= 5 and 0 <= x2 <= 9 with accuracy correct to
(the default) 3 decimal digits.
In case f1 and f2 are real and imaginary parts of a single
function f, the result will be all the (real and complex) roots of the
function f(x) in the specified rectangle.
For example,
comroots[{Exp[2x1] Cos[2x2] + 3,
Exp[2x1] Sin[2x2]}, {x1, 0, 1}, {x2, 0, 30}]
finds all the (complex) roots of Tanh(x) = 2 on the rectangle
[0, 1] x [0, 30].
In case f(x) is a POLYNOMIAL with real or complex coefficients, after
running the accompanying pacakge reim. m (from Maeder),
cxroots[f(x), {x1, a, b}, {x2, c, d}, tol]
will give all the real and complex roots of the polynomial f in the
rectangle [a, b] x [c, d].
For example, try
cxroots[x^3 +1, {x1, -2, 1}, {x2, -1, 1}]
(Sorry for the repetitions in the output...).
The output will be displayed in terms of the Mathematica built-in
function RealInterval. {RealInterval[{x1,x2}], RealInterval[{y1,y2}]
will stand for the root which is contained in the closed
rectangle [x1, x2] x [y1, y2]. This output format is preferable
since Mathematica can do arithmetic with RealInterval objects.
This is a rather slow programme. Faster versions, using
Secant and Newton's methods, are also available.
The built-in function RealInterval is not yet bug-free.
I would be interested to know if you come across any anomolies
related to this package.
Yilmaz Akyildiz
akyildiz at maths.oxford.ac.uk
*)
BeginPackage["ReIm`"];
Begin["`Private`"];
protected = Unprotect[Re, Im]
ri[x_,y_] := RealInterval[{x, y}]
Re[x_] := x /; Im[x] == 0
Re[x_ + y_] := Re[x] + Re[y]
Im[x_ + y_] := Im[x] + Im[y]
Re[x_ y_] := Re[x] Re[y] - Im[x] Im[y]
Im[x_ y_] := Re[x] Im[y] + Im[x] Re[y]
reim[f_] := {Re[Expand[f/.x->x+I y]], Im[Expand[f/.x->x+I y]]} /.
{Im[x]->0,Im[y]->0,Re[x]->x,Re[y]->y,
Im[x^n_]->0,Im[y^n_]->0,Re[x^n_]->x^n,Re[y^n_]->y^n}
Protect[ Evaluate[protected] ]
EndPackage[ ]
ri[x_, y_] := RealInterval[{x, y}]
leftup[{ri[x1_, x2_], ri[y1_, y2_]}] :=
{ri[x1, (x1+x2)/2.], ri[(y1+y2)/2., y2]} // N
leftdown[{ri[x1_, x2_], ri[y1_, y2_]}] :=
{ri[x1, (x1+x2)/2.], ri[y1, (y1+y2)/2.]} // N
rightup[{ri[x1_, x2_], ri[y1_, y2_]}] :=
{ri[(x1+x2)/2., x2], ri[(y1+y2)/2., y2]} // N
rightdown[{ri[x1_, x2_], ri[y1_, y2_]}] :=
{ri[(x1+x2)/2., x2], ri[y1, (y1+y2)/2.]} // N
quadrasect[rect_] := {leftup[rect], leftdown[rect],
rightup[rect], rightdown[rect]}
containszero[{ri[x1_, x2_], ri[y1_, y2_]}] :=
(x1 <= 0 <= x2) && (y1 <= 0 <= y2)
F[{f1_,f2_}, {a_RealInterval,b_RealInterval}] :=
{f1, f2} /. {x1->a, x2->b}
comroots[{f1_, f2_}, {x1_, a_, b_}, {x2_, c_, d_}, tol_:0.0001] :=
With[{step = Select[Flatten[(quadrasect/@#),1],
containszero[N[F[{f1,f2},#]]]&]& },
Nest[step, N[{{ri[a, b], ri[c, d]}}],
Ceiling[Log[2, N[Max[b-a,d-c]/tol]]]]
]
cxroots[f_, {x1_, a_, b_}, {x2_, c_, d_}, tol_:0.00001] :=
comroots[reim[f]/.{x->x1, y->x2},
{x1, a, b}, {x2, c, d}, tol]