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