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

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] :=
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]

```

• Prev by Date: Re: HypergeometricPFQ
• Next by Date: Math typesetting in plot labels
• Previous by thread: Save and big arrays