(************** Content-type: application/mathematica ************** CreatedBy='Mathematica 4.2' Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 15463, 699]*) (*NotebookOutlinePosition[ 16114, 722]*) (* CellTagsIndexPosition[ 16070, 718]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["LaguerreNRoot", "Title"], Cell[TextData[StyleBox["T. M. Jonassen & R. Ghulghazaryan", FontSlant->"Italic"]], "Text", FontSize->14], Cell[TextData[{ "The following example implements the ", StyleBox["laguer()", FontWeight->"Bold"], " routine from NRIC as an engine for a rootfinder for polynomials in one \ complex variable with complex coefficients. The function, called ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], ", serves the same purpose as ", StyleBox["FindRoot", FontWeight->"Bold"], ". ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], " is limited to polynomials while ", StyleBox["FindRoot", FontWeight->"Bold"], " is not, it works for arbitrarly equations. The syntax is also slightly \ different. To find a root for the polynomial ", StyleBox["poly==0", FontSlant->"Italic"], ", ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], " uses the syntax ", StyleBox["LaguerreNRoot[", FontWeight->"Bold"], StyleBox["poly", FontSlant->"Italic"], ",{", StyleBox["var, guess", FontSlant->"Italic"], "}", StyleBox["]", FontWeight->"Bold"], " while ", StyleBox["FindRoot", FontWeight->"Bold"], " uses ", StyleBox["FindRoot[", FontWeight->"Bold"], StyleBox["poly==0", FontSlant->"Italic"], ",{", StyleBox["var, guess", FontSlant->"Italic"], "}]. The initial guess in ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], " may be real or complex." }], "Text"], Cell[CellGroupData[{ Cell["Utillity functions", "Section"], Cell[TextData[{ "We will need a few simple functions in order to define LaguerreNRoot. The \ function of these are easily understood and not commented here. They are same \ as used for ", StyleBox["LaguerreNSolve", FontWeight->"Bold"], ", in fact we need one less:" }], "Text"], Cell[BoxData[{ \(\(rtc[l_] := {Re[l], Im[l]};\)\), "\[IndentingNewLine]", \(\(makelist[poly_, var_] := Flatten[Map[rtc, CoefficientList[poly, var]]];\)\), "\[IndentingNewLine]", \(\(tocomp[{\[Alpha]_, \[Beta]_}] := \[Alpha] + \[ImaginaryI]\ \[Beta];\)\ \)}], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["External function", "Section"], Cell[TextData[{ "We must install the external program ", StyleBox["oneroot", FontWeight->"Bold"], ". Note that this program must be running under use of ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], "." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(link = Install["\"]\)], "Input"], Cell[OutputFormData["\<\ LinkObject[\"/home/torejo/mathbin/oneroot\", 2, 2]\ \ \>", "\<\ LinkObject[/home/torejo/mathbin/oneroot, 2, 2]\ \>"], "Output"] }, Open ]], Cell["\<\ The following terminates the external program (do not use this now \ -:)):\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Uninstall[link]\)], "Input"], Cell[OutputFormData["\<\ \"/home/torejo/mathbin/oneroot\"\ \>", \ "\<\ /home/torejo/mathbin/oneroot\ \>"], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["LaguerreNRoot definition", "Section"], Cell[TextData[{ "The construction of ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], " is straight forward. The main function is the external function ", StyleBox["getoneroot", FontWeight->"Bold"], " defined by the external program ", StyleBox["oneroot", FontWeight->"Bold"], ". The rest of the code is simply formatting stuff. The output is a rule, \ the same as in ", StyleBox["FindRoot", FontWeight->"Bold"], "." }], "Text"], Cell[BoxData[ \(LaguerreNRoot[ poly_, {var_, guess_}] := {Rule[var, Chop[tocomp[ Evaluate[ getoneroot[Evaluate[N[makelist[poly, var]]], N[rtc[guess]]]]]]]}\)], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Limitations", "Section"], Cell[TextData[{ "The polynomial must have numerical coefficients, and be a polynomial in \ one variable only. The polynomial must have non-negative integer exponents. \ The code has not been extensively tested, and ", StyleBox["the error handeling could be better", FontSlant->"Italic"], "." }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Examples", "Section"], Cell[TextData[{ "Here are some simple use of ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], ". As in ", StyleBox["FindRoot", FontWeight->"Bold"], ", the solution is returned as a rule depending on the initial guess:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[x\^2 + x - 2, {x, 1}]\)], "Input"], Cell[OutputFormData["\<\ {x -> 1.}\ \>", "\<\ {x -> 1.}\ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[x\^2 + x - 2, {x, \(-1\)}]\)], "Input"], Cell[OutputFormData["\<\ {x -> -2.}\ \>", "\<\ {x -> -2.}\ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[x\^7 - 4 x\^3 + 2, {x, 3}]\)], "Input"], Cell[OutputFormData["\<\ {x -> 1.3327731271582208}\ \>", "\<\ {x -> \ 1.33277}\ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[x\^7 - 4 x\^3 + 2, {x, 1}]\)], "Input"], Cell[OutputFormData["\<\ {x -> 0.827357153301707}\ \>", "\<\ {x -> \ 0.827357}\ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[x\^7 - 4 x\^3 + 2, {x, 0}]\)], "Input"], Cell[OutputFormData["\<\ {x -> -0.4055983097396118 - \ 0.6655184609220427*I}\ \>", "\<\ {x -> -0.405598 - 0.665518 I}\ \>"], "Output"] }, Open ]], Cell["\<\ The argument can be anything that evaluates to a polynomial of the \ desired form:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[ Sum[Random[Integer, {\(-2\), 2}]\ \[Zeta]\^i, {i, 0, 10}], {\[Zeta], 3}]\)], "Input"], Cell[OutputFormData["\<\ {\[Zeta] -> -0.1988271030065306 + \ 0.9020591278825711*I}\ \>", "\<\ {\[Zeta] -> -0.198827 + 0.902059 I}\ \>"], \ "Output"] }, Open ]], Cell[TextData[{ StyleBox["LaguerreNRoot", FontWeight->"Bold"], " fails if the polynomial contains non-numerical coefficients, or if the \ guess is not a real or complex number:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[s\^4 + \[Alpha]\ s\ - 1, {s, \[ImaginaryI]}]\)], "Input"], Cell[OutputFormData["\<\ {s -> tocomp[$Failed]}\ \>", "\<\ {s -> \ tocomp[$Failed]}\ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[s\^4 + \ s\ - 1, {s, {0, 1}}]\)], "Input"], Cell[OutputFormData["\<\ {s -> tocomp[$Failed]}\ \>", "\<\ {s -> \ tocomp[$Failed]}\ \>"], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[TextData[{ StyleBox["FindRoot", FontWeight->"Plain"], " comparison" }], "Section"], Cell[TextData[{ "In simple cases ", StyleBox["LaguerreNRoot", FontWeight->"Bold"], " and ", StyleBox["FindRoot", FontWeight->"Bold"], " returns the same answer:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[ x\^7 - 4 x\^3 + 2 \[ImaginaryI]\ x\^2 - x + 1, {x, 0}]\)], "Input"], Cell[OutputFormData["\<\ {x -> 0.46280204108650763 + \ 0.1157730755483299*I}\ \>", "\<\ {x -> 0.462802 + 0.115773 I}\ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(FindRoot[ x\^7 - 4 x\^3 + 2 \[ImaginaryI]\ x\^2 - x + 1 \[Equal] 0, {x, 0}]\)], "Input"], Cell[OutputFormData["\<\ {x -> 0.46280202044013335 + 0.11577307663251217*I}\ \ \>", "\<\ {x -> 0.462802 + 0.115773 I}\ \>"], "Output"] }, Open ]], Cell["In other cases we get different answers:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[Nest[\((#\^2 - 5)\) &, x, 5], {x, 1}]\)], "Input"], Cell[OutputFormData["\<\ {x -> 1.3453038108409354}\ \>", "\<\ {x -> 1.3453}\ \ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(FindRoot[Nest[\((#\^2 - 5)\) &, x, 5] \[Equal] 0, {x, 1}]\)], "Input"], Cell["\<\ FindRoot::cvnwt: Newton's method failed to converge to the prescribed accuracy after 15 iterations.\ \>", "Message"], Cell[OutputFormData["\<\ {x -> 1.4845340746948277}\ \>", "\<\ {x -> \ 1.48453}\ \>"], "Output"] }, Open ]], Cell["\<\ It seems that LaguerreNRoot does not compute correctly, the \ following is wrong:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[Nest[\((#\^2 - 5)\) &, x, 8], {x, 11}]\)], "Input"], Cell[OutputFormData["\<\ {x -> 8.542198086440466}\ \>", "\<\ {x -> \ 8.5422}\ \>"], "Output"] }, Open ]], Cell["Changing guess it does not compute at all:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LaguerreNRoot[Nest[\((#\^2 - 5)\) &, x, 10], {x, 1}]\)], "Input"], Cell[OutputFormData["\<\ {x -> tocomp[$Failed]}\ \>", "\<\ {x -> \ tocomp[$Failed]}\ \>"], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(FindRoot[Nest[\((#\^2 - 5)\) &, x, 10] \[Equal] 0, {x, 1}]\)], "Input"], Cell["\<\ FindRoot::cvnwt: Newton's method failed to converge to the prescribed accuracy after 15 iterations.\ \>", "Message"], Cell[OutputFormData["\<\ {x -> 1.0372094741728022}\ \>", "\<\ {x -> \ 1.03721}\ \>"], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["References", "Section"], Cell[TextData[{ StyleBox["NRIC W. H. Press, S. A. Teukolsky, W. T. Vetterling & B. P. \ Flannery: ", FontSize->10], StyleBox["Numerical Recipes in C", FontSize->10, FontSlant->"Italic"], StyleBox[". Cambridge Univerity Press. ", FontSize->10], StyleBox["1992", FontSize->10, FontWeight->"Bold"] }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Compilation", "Section"], Cell[TextData[{ "On a Linux platform the ", StyleBox["MathLink", FontSlant->"Italic"], " program is compiled by:" }], "Text"], Cell["mcc -o oneroot -O2 -march=i686 oneroot.tm oneroot.c -lm -lML", "Program"], Cell[TextData[{ "My suggestion is to keep all working ", StyleBox["MathLink", FontSlant->"Italic"], " programs in one directory. On my machine I have used the directory ", StyleBox["/home/torejo/mathbin", FontFamily->"Lucidatypewriter"], " for private programs." }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["MathLink template file", "Section"], Cell[TextData[{ "This is the ", StyleBox["MathLink", FontSlant->"Italic"], " template file, here called ", StyleBox["oneroot.tm", FontWeight->"Bold", FontSlant->"Italic"], ":" }], "Text"], Cell["\<\ :Begin: :Function: getoneroot :Pattern: getoneroot[c_List,g_List] :Arguments: {c,g} :ArgumentTypes: {RealList,RealList} :ReturnType: Manual :End:\ \>", "Program"] }, Open ]], Cell[CellGroupData[{ Cell["C code for external functions", "Section"], Cell[TextData[{ "This is the C code, here called ", StyleBox["oneroot.c", FontWeight->"Bold", FontSlant->"Italic"], ":" }], "Text"], Cell["\<\ /* MathLink program for LaguerreNRoot */ /* Compile with : mcc -o oneroot -O2 -march=i686 oneroot.tm oneroot.c -lm \ -lML */ #include #include #include \"/usr/local/Wolfram/Mathematica/4.2/AddOns/MathLink/DeveloperKit/\ Linux/CompilerAdditions/mathlink.h\" #include #include /* #define's for laguerre() */ #define EPSS 1.0e-7 #define MR 8 #define MT 10 #define MAXIT (MT*MR) /* Macro used in laguerre() */ static double maxarg1, maxarg2; #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2) ? (maxarg1) : \ (maxarg2)) /* Find one complex root, a root finder from NRIC */ void laguerre(complex *a, int m, complex *x, int *its) { \tint iter,j; \tdouble abx,abp,abm,err; \tcomplex dx,x1,b,d,f,g,h,sq,gp,gm,g2; \tstatic double frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0}; \tfor (iter=1;iter<=MAXIT;iter++) { \t\t*its=iter; \t\tb=a[m]; \t\terr=cabs(b); \t\td=0.0; \t\tf=0.0; \t\tabx=cabs(*x); \t\tfor (j=m-1;j>=0;j--) { \t\t\tf=(*x)*f+d; \t\t\td=(*x)*d+b; \t\t\tb=(*x)*b+a[j]; \t\t\terr=cabs(b)+abx*err; \t\t} \t\terr *= EPSS; \t\tif (cabs(b) <= err) return; \t\tg=d/b; \t\tg2=g*g; \t\th=g2-(2.0*(f/b)); \t\tsq=csqrt((double) (m-1)*((double) m*h-g2)); \t\tgp=g+sq; \t\tgm=g-sq; \t\tabp=cabs(gp); \t\tabm=cabs(gm); \t\tif (abp < abm) gp=gm; \t\tdx = FMAX(abp,abm) > 0.0 ? ((double)m+0.0*I)/gp \t\t\t: exp(log(1+abx))*(cos((double)iter)+I*sin((double)iter)); \t\tx1=(*x)-dx; \t\tif (creal(*x) == creal(x1) && cimag(*x) == cimag(x1)) return; \t\tif (iter % MT) *x=x1; \t\telse *x = *x-(frac[iter/MT]-dx); \t} \tfprintf(stderr,\"too many iterations in laguerre\\n\"); \treturn; } void getoneroot(double *c, long length, double *guess, long lg) { complex *a; complex *x; int *iter; double *retarr; long dp; int i,val; int deg; if(lg != 2) { fprintf(stderr,\"Initial guess is not Real or Complex\\n\"); exit(1); } deg = (int)(length/2) - 1; dp = 2; /* Allocate space for \"a\", \"x\" and \"retarr\" */ a = (complex *)malloc((deg+1)*sizeof(complex)); x = (complex *)malloc(sizeof(complex)); retarr = (double *)malloc(2*sizeof(double)); iter = (int *)malloc(sizeof(int)); /* Fill the array a with correct values */ for(i=0; i", "Program"] }, Open ]] }, Open ]] }, FrontEndVersion->"4.2 for X", ScreenRectangle->{{0, 1600}, {0, 1200}}, WindowSize->{1497, 579}, WindowMargins->{{Automatic, 21}, {Automatic, 89}}, Magnification->1.5 ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1776, 53, 30, 0, 169, "Title"], Cell[1809, 55, 109, 2, 50, "Text"], Cell[1921, 59, 1353, 50, 121, "Text"], Cell[CellGroupData[{ Cell[3299, 113, 37, 0, 87, "Section"], Cell[3339, 115, 287, 7, 71, "Text"], Cell[3629, 124, 303, 6, 85, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[3969, 135, 36, 0, 87, "Section"], Cell[4008, 137, 236, 8, 46, "Text"], Cell[CellGroupData[{ Cell[4269, 149, 83, 1, 39, "Input"], Cell[4355, 152, 152, 5, 36, "Output"] }, Open ]], Cell[4522, 160, 98, 3, 46, "Text"], Cell[CellGroupData[{ Cell[4645, 167, 48, 1, 39, "Input"], Cell[4696, 170, 116, 5, 36, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[4861, 181, 43, 0, 87, "Section"], Cell[4907, 183, 463, 15, 71, "Text"], Cell[5373, 200, 240, 6, 62, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[5650, 211, 30, 0, 87, "Section"], Cell[5683, 213, 314, 7, 71, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[6034, 225, 27, 0, 87, "Section"], Cell[6064, 227, 249, 8, 46, "Text"], Cell[CellGroupData[{ Cell[6338, 239, 68, 1, 47, "Input"], Cell[6409, 242, 72, 4, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[6518, 251, 73, 1, 47, "Input"], Cell[6594, 254, 74, 4, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[6705, 263, 74, 1, 47, "Input"], Cell[6782, 266, 95, 5, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[6914, 276, 74, 1, 47, "Input"], Cell[6991, 279, 95, 5, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[7123, 289, 74, 1, 47, "Input"], Cell[7200, 292, 134, 5, 36, "Output"] }, Open ]], Cell[7349, 300, 106, 3, 46, "Text"], Cell[CellGroupData[{ Cell[7480, 307, 135, 3, 47, "Input"], Cell[7618, 312, 148, 6, 36, "Output"] }, Open ]], Cell[7781, 321, 197, 5, 46, "Text"], Cell[CellGroupData[{ Cell[8003, 330, 92, 1, 47, "Input"], Cell[8098, 333, 100, 5, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[8235, 343, 77, 1, 47, "Input"], Cell[8315, 346, 100, 5, 36, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[8464, 357, 95, 4, 87, "Section"], Cell[8562, 363, 191, 8, 46, "Text"], Cell[CellGroupData[{ Cell[8778, 375, 110, 2, 47, "Input"], Cell[8891, 379, 133, 5, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[9061, 389, 125, 3, 47, "Input"], Cell[9189, 394, 134, 5, 36, "Output"] }, Open ]], Cell[9338, 402, 56, 0, 46, "Text"], Cell[CellGroupData[{ Cell[9419, 406, 84, 1, 47, "Input"], Cell[9506, 409, 94, 5, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[9637, 419, 90, 1, 47, "Input"], Cell[9730, 422, 135, 4, 78, "Message"], Cell[9868, 428, 95, 5, 36, "Output"] }, Open ]], Cell[9978, 436, 105, 3, 46, "Text"], Cell[CellGroupData[{ Cell[10108, 443, 85, 1, 47, "Input"], Cell[10196, 446, 93, 5, 36, "Output"] }, Open ]], Cell[10304, 454, 58, 0, 46, "Text"], Cell[CellGroupData[{ Cell[10387, 458, 85, 1, 47, "Input"], Cell[10475, 461, 100, 5, 36, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[10612, 471, 91, 1, 47, "Input"], Cell[10706, 474, 135, 4, 78, "Message"], Cell[10844, 480, 95, 5, 36, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[10988, 491, 29, 0, 87, "Section"], Cell[11020, 493, 341, 12, 44, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[11398, 510, 30, 0, 87, "Section"], Cell[11431, 512, 135, 5, 46, "Text"], Cell[11569, 519, 79, 0, 57, "Program"], Cell[11651, 521, 293, 8, 47, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[11981, 534, 41, 0, 87, "Section"], Cell[12025, 536, 209, 9, 46, "Text"], Cell[12237, 547, 172, 8, 177, "Program"] }, Open ]], Cell[CellGroupData[{ Cell[12446, 560, 48, 0, 87, "Section"], Cell[12497, 562, 146, 6, 46, "Text"], Cell[12646, 570, 2789, 125, 2457, "Program"] }, Open ]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)