Re: ContourFindRoot
- To: mathgroup at smc.vnet.net
- Subject: [mg71327] Re: ContourFindRoot
- From: "dimitris" <dimmechan at yahoo.com>
- Date: Wed, 15 Nov 2006 06:43:23 -0500 (EST)
- References: <ejc4rs$6nr$1@smc.vnet.net>
In the previous post I forrgot to mention any kind of references. So here they are http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/f042e78929ddb078/0a9c893ced104f6f?lnk=st&q=FindRoot&rnum=3#0a9c893ced104f6f http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/2587d839a40223bb/1a2c1c290d3c4b26?lnk=gst&q=FindRoot+David+Park&rnum=1#1a2c1c290d3c4b26 http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/9836d5644ef5b0a1/49569b37ba055a99?lnk=gst&q=FindRoot&rnum=7#49569b37ba055a99 http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/7612c1825cd938da/f5c53fdd2c7b737e?lnk=gst&q=Find+crossings&rnum=3#f5c53fdd2c7b737e http://www.amazon.com/Mathematica-Action-Stan-Wagon/dp/0387986847/sr=1-2/qid=1162830266/ref=sr_1_2/103-5414680-6091022?ie=UTF8&s=books http://www.amazon.com/Complex-Analysis-MATHEMATICA-William-Shaw/dp/0521836263/sr=8-1/qid=1163459323/ref=pd_bbs_sr_1/103-5414680-6091022?ie=UTF8&s=books Dimitris dimitris wrote: > 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