Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

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

Search the Archive

ContourFindRoot

  • To: mathgroup at smc.vnet.net
  • Subject: [mg71302] ContourFindRoot
  • From: "dimitris" <dimmechan at yahoo.com>
  • Date: Tue, 14 Nov 2006 05:06:19 -0500 (EST)

Below we construct a routine, which will take two equations {f[x,y]=0,
g[x,y]=0} and return all the simultaneous solutions (only the
crossings; not tangencies) in a rectangle.

The basic steps of the routine are the following:

First, we generate a contour plot of f[x,y]=0, go into the graphics
objects that results, and gather up the data the data making up the
curves. Then, on each curve, we evaluate g. The places where g changes
sign correspond to solutions; perhaps not very accurate. We next send
them to FindRoot in order to get exact solution.

----------------------------------------------------------------------------------------------------------------------

Clear["Global`*"]

Off[General::spell1]

CurveFindRoot::usage =
"CurveFindRoot[{f,g},{x,xmin,xmax},{y,ymin,ymax},curve\
,points,opts]\nfinds all nontangential solutions to {f=0, g=0}, on a
contour \
line.";

ContourFindRoot::usage =
"ContourFindRoot[{f,g},{x,xmin,xmax},{y,ymin,ymax},p\
oints,opts]\nfinds all nontangential solutions to {f=0, g=0}. in a
rectangle \
in the 2D domain";

ContourPlotFunctionRoots::usage =
"ContourPlotFunctionRoots[{f,g},{x,xmin,xma\
x},{y,ymin,ymax},points,roots,optsplot]\ngives the graph of curves and
\
zeros.";

CurveFindRoot[{f_, g_}, {x_, x0_, x1_}, {y_, y0_, y1_}, curve_,
points_,
   opts___] := Block[{contourdata, pos, seeds, gg},
   contourdata = Cases[Graphics[ContourPlot[f, {x, x0, x1}, {y, y0,
y1},
        Contours -> {0}, ContourShading -> False, PlotPoints -> points,

        DisplayFunction -> Identity]], Line[a_] :> a, Infinity];
    gg = Function[{x, y}, Evaluate[g]];
    pos = Position[Partition[Apply[gg, contourdata[[curve]], {1}], 2,
1],
      _?(Times @@ #1 < 0 & ), {1}]; seeds = contourdata[[curve,
      Flatten[pos]]]; (FindRoot[{f == 0, g == 0}, {x, #1[[1]]}, {y,
#1[[2]]},
       opts] & ) /@ seeds]

ContourFindRoot[{f_, g_}, {x_, x0_, x1_}, {y_, y0_, y1_}, points_,
   opts___] := Block[{lengthcontourdata, roots},
   lengthcontourdata = Length[Cases[Graphics[ContourPlot[f, {x, x0,
x1},
         {y, y0, y1}, Contours -> {0}, ContourShading -> False,
         PlotPoints -> points, DisplayFunction -> Identity]], Line[a_]
:> a,
       Infinity]]; roots = (CurveFindRoot[{f, g}, {x, x0, x1}, {y, y0,
y1},
        #1, points, opts] & ) /@ Range[lengthcontourdata];
    {x, y} /. Flatten[DeleteCases[roots, {}], 1]]

ContourPlotFunctionRoots[{f_, g_}, {x_, x0_, x1_}, {y_, y0_, y1_},
points_,
   roots_, optsplot___] :=
  Show[MapThread[ContourPlot[#1, {x, x0, x1}, {y, y0, y1}, optsplot,
      PlotPoints -> points, Contours -> {0}, ContourShading -> False,
      DisplayFunction -> Identity, ContourStyle -> {Hue[#2/3]}] & ,
    {{f, g}, {1, 2}}], Graphics[
    ({Red, AbsolutePointSize[6], Point[#1]} & ) /@ roots],
   DisplayFunction -> $DisplayFunction]

---------------------------------------------------------------------------------------------------------------------

Here are some applications on finding the roots of a function in a
rectangle
in the complex domain.

Clear[f]

f[z_] := Sin[z^3] + Cos[z^2] - z

ref[x_, y_] = Simplify[Re[ComplexExpand[f[z] /. z -> x + I*y]],
   Element[{x,y},Reals]]
imf[x_, y_] = Simplify[Im[ComplexExpand[f[z] /. z -> x + I*y]],
   Element[{x,y},Reals]]
-x + Cos[x^2 - y^2]*Cosh[2*x*y] + Cosh[3*x^2*y - y^3]*Sin[x^3 -
3*x*y^2]
-y - Sin[x^2 - y^2]*Sinh[2*x*y] + Cos[x^3 - 3*x*y^2]*Sinh[3*x^2*y -
y^3]

sols = ContourFindRoot[{ref[x, y], imf[x, y]}, {x, -2, 2}, {y, -2, 2},
100]

ContourPlotFunctionRoots[{ref[x, y], imf[x, y]}, {x, -2, 2}, {y, -2,
2}, 100, sols];

Here is another example adopted from a private communication with David
Park

Clear[f]

f[z_] := Sin[Cos[z]] - z

ref[x_, y_] = Simplify[Re[ComplexExpand[f[z] /. z -> x + I*y]],
  Element[{x,y},Reals]]
imf[x_, y_] = Simplify[Im[ComplexExpand[f[z] /. z -> x + I*y]],
   Element[{x,y},Reals]]
-x + Cosh[Sin[x]*Sinh[y]]*Sin[Cos[x]*Cosh[y]]
-y - Cos[Cos[x]*Cosh[y]]*Sinh[Sin[x]*Sinh[y]]

Timing[sols = ContourFindRoot[{ref[x, y], imf[x, y]}, {x, -Pi, Pi},
    {y, -Pi, Pi}, 200]]

ContourPlotFunctionRoots[{ref[x, y], imf[x, y]}, {x, -Pi, Pi}, {y, -Pi,
Pi},
   200, sols];

Timing[sols = ContourFindRoot[{ref[x, y], imf[x, y]}, {x, -Pi, Pi},
    {y, Pi, 5*(Pi/3)}, 200]]

ContourPlotFunctionRoots[{ref[x, y], imf[x, y]}, {
    x, -Pi, Pi}, {y, Pi, 5Pi/3}, 200, sols, ImageSize -> {600, 400}];

Here is a system of two equations f[x,y]==0&&g[x,y]==0

f[x_, y_] := Cosh[x + y] + x^2*y^2 - 7*y*x^3
g[x_, y_] := Sinh[x + y] + x*y - 2*x

Timing[sols = ContourFindRoot[{f[x, y], g[x, y]}, {x, -10, 10}, {y,
-10, 10}, 200]]

ContourPlotFunctionRoots[{f[x, y], g[x, y]}, {x, -10, 10}, {y, -10,
10}, 200, sols]

I would appreciate a lot any kinds of comments for my code, possible
modifications,
suggestions for extension to treat more cases etc.

I would further appreciate if you experiment with it with functions of
your own
and write down the results.

Best Regards
Dimitris


  • Prev by Date: Re: Re: finding the (v,w) weighted degree of a polynomial
  • Next by Date: Scratching an image!
  • Previous by thread: Re: FindUnknowns[ ] (Mathematica wish list)
  • Next by thread: Re: ContourFindRoot