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

Re: 2 dimension Newton Raphson

  • To: mathgroup at smc.vnet.net
  • Subject: [mg71276] Re: 2 dimension Newton Raphson
  • From: "dimitris" <dimmechan at yahoo.com>
  • Date: Sun, 12 Nov 2006 06:49:01 -0500 (EST)
  • References: <ej1q7l$e2l$1@smc.vnet.net>

Here is the solutions using Solve and NSolve

eq1 = (x - 4)^2 + (y - 4)^2 == 5;
eq2 = x^2 + y^2 == 16;

Solve[{eq1, eq2}, {x, y}]
{{x -> (1/16)*(43 - Sqrt[199]), y -> (1/16)*(43 + Sqrt[199])}, {x ->
(1/16)*(43 + Sqrt[199]), y -> (1/16)*(43 - Sqrt[199])}}

N[%]
{{x -> 1.8058290012708822, y -> 3.569170998729118}, {x ->
3.569170998729118, y -> 1.8058290012708822}}

NSolve[{eq1, eq2}, {x, y}]
{{x -> 1.8058290012708826, y -> 3.5691709987291174}, {x ->
3.569170998729118, y -> 1.8058290012708822}}

Here is another approach using my rootine ContourFindRoot.

I constructed a routine, which  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 send them to
FindRoot in order to get exact
solution.

I used the folloing references:

http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/9836d5644ef5b0a1/49569b37ba055a99?lnk=gst&q=FindRoot&rnum=7#49569b37ba055a99

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://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/2587d839a40223bb/1a2c1c290d3c4b26?lnk=gst&q=Park+findRoot&rnum=1#1a2c1c290d3c4b26

http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/7612c1825cd938da/f5c53fdd2c7b737e?lnk=gst&q=Find+crossings&rnum=3#f5c53fdd2c7b737e

Here is the code:

Clear["Global`*"]

Off[General::spell1]

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

ContourFindRoot::usage =
   "ContourFindRoot[{f,g},{x,xmin,xmax},{y,ymin,ymax},points,opts]
finds all nontangential solutions to \n\t\t{f=0, g=0}.";

ContourPlotFunctionRoots::usage =

"ContourPlotFunctionRoots[{f,g},{x,xmin,xmax},{y,ymin,ymax},roots,optsplot]
gives 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 is an application to your problem

f[x_, y_] := (x - 4)^2 + (y - 4)^2 - 5
g[x_, y_] := x^2 + y^2 - 16

sols = ContourFindRoot[{f[x, y], g[x, y]}, {x, -5, 5}, {y, -5, 5}, 200]
{{3.569170998729118, 1.8058290012708822}, {1.8058290012708818,
3.5691709987291182}}

or if you wish in rules format

({x -> #1[[1]], y -> #1[[2]]} & ) /@ ContourFindRoot[{f[x, y], g[x,
y]}, {x, -5, 5}, {y, -5, 5}, 200]
{{x -> 3.569170998729118, y -> 1.8058290012708822}, {x ->
1.8058290012708818, y -> 3.5691709987291182}}

And here is a contour plot of (x - 4)^2 + (y - 4)^2 - 5=0 and x^2 + y^2
- 16=0 along the
roots of their system

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

Regards
Dimitris Anagnostou


ms z wrote:
> I have tried to solve the roots of the simultaneous nonlinear equations
> (x-4)^2 + (y-4)^2 = 5
> x^2 + y^2 = 16
>
> by writing this function:
>
> nr2method[xl1_, xl2_, es1_] :=
>   Block[{x1, x2, ea, es, x1new, u, v},
>     u = (x1 - 4)^2 + (x2 - 4)^2 - 5;
>     v = x1^2 + x2^2 - 16;
>     ea = 100; es = es1;
>     For[i = 1, ea > es, i++,
>       (x1new[x1_, x2_] = x1 - (u*D[
>     v, x2] - v*D[u, x2])/(D[u, x1]*D[v, x2] - D[u, x2]*D[v, x1]);
>         If[i == 1, x1 = xl1, x1 = b];
>         x2 = xl2;
>         b = x1new[x1, x2];
>         ea = Abs[(b - x1)/b 100];
>         Clear[x1, x2, x1new];)];
>     ea = 100; es = es1;
>     For[i = 1, ea > es, i++,
>       (x2new[x1_, x2_] = x2 - (v*D[u, x1] - u*D[v, x1])/(D[u, x1]*D[v, x2] -
>         D[u, x2]*D[v, x1]);
>         If[i == 1, x2 = xl2, x2 = c];
>         x1 = xl1;
>         c = x2new[x1, x2];
>         ea = Abs[(c - x2)/c 100];
>         Clear[x1, x2, x2new];)];
>     Print["The value of x1 is ", b];
>     Print["The value of x2 is ", c];]
>
> Is this function a good one? Is there a way to make this function simpler?
>


  • Prev by Date: Re: are there any methods of figuring out how "large" a piece of typeset textual data will be?
  • Next by Date: Re: Function to solve polynomial
  • Previous by thread: RE: 2 dimension Newton Raphson
  • Next by thread: Re: 2 dimension Newton Raphson