{VERSION 3 0 "IBM INTEL NT" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 1 10 128 0 0 1 0 2 0 0 1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 1 11 0 0 0 1 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Comment" 2 18 "" 0 0 0 128 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "2D Input" 2 19 "" 0 0 128 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "2 D Output" 2 20 "" 0 0 0 0 128 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "Non-propor tional" -1 256 "Courier" 1 11 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 0 2 0 0 0 0 0 0 }0 1 0 -1 6 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 24 6 0 0 0 0 0 0 -1 0 }{PSTYLE " Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 15 0 0 0 0 0 0 0 0 0 0 0 0 0 } 0 0 0 -1 6 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE " " -1 -1 "" 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 12 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Bullet Item" 0 15 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 3 3 0 0 0 0 0 0 15 2 }{PSTYLE "Title " 0 18 1 {CSTYLE "" -1 -1 "Times" 1 22 0 128 128 1 0 1 0 0 0 0 0 0 0 } 3 1 0 -1 14 12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "Times" 1 16 0 0 128 1 0 1 0 0 0 0 0 0 0 }3 1 0 -1 8 0 0 0 0 0 0 0 256 0 }{PSTYLE "Address" 0 256 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 128 1 1 2 0 0 0 0 0 0 0 }3 1 0 -1 6 -1 0 0 0 0 0 0 258 0 }{PSTYLE "D ate" 0 257 1 {CSTYLE "" -1 -1 "" 1 14 0 0 128 1 2 0 0 0 0 0 0 0 0 }3 0 0 -1 8 6 0 0 0 0 0 0 -1 0 }{PSTYLE "Address - email" -1 258 1 {CSTYLE "" -1 -1 "Verdana" 1 10 0 0 128 1 2 0 0 0 0 0 0 0 0 }3 0 0 -1 4 -1 0 0 0 0 0 0 257 0 }} {SECT 0 {PARA 18 "" 0 "" {TEXT -1 49 "A Symbolic Newton-Raphson Method of Finding Roots" }}{EXCHG {PARA 19 "" 0 "" {TEXT -1 15 "Marc A. Muri son" }}{PARA 256 "" 0 "" {TEXT -1 59 "Astronomical Applications Depart ment\nU.S. Naval Observatory" }}{PARA 258 "" 0 "" {TEXT -1 57 "murison @aa.usno.navy.mil\nhttp://aa.usno.navy.mil/murison/" }}{PARA 257 "" 0 "" {TEXT -1 13 "December 1996" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 15 "1. Introduction" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 " Suppose we have a function " }{XPPEDIT 18 0 "f(x)" "6#-%\"fG6#%\"xG" }{TEXT -1 62 " for which we'd like to find a root. At some arbitrary point " } {XPPEDIT 18 0 "x=a" "6#/%\"xG%\"aG" }{TEXT -1 38 ", the equation of th e line tangent to " }{XPPEDIT 18 0 "f(x)" "6#-%\"fG6#%\"xG" }{TEXT -1 4 " is " }}{PARA 0 "" 0 "" {TEXT -1 7 "(1) " }{XPPEDIT 18 0 "y = di ff(f(a),x)*x + b" "6#/%\"yG,&*&-%%diffG6$-%\"fG6#%\"aG%\"xG\"\"\"F.F/F /%\"bGF/" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 6 "where " } {TEXT 257 1 "b" }{TEXT -1 8 " is the " }{TEXT 258 1 "y" }{TEXT -1 53 " intercept and we use the notational convention that " }{XPPEDIT 18 0 "diff(f(a),x)" "6#-%%diffG6$-%\"fG6#%\"aG%\"xG" }{TEXT -1 22 " is the \+ derivative of " }{XPPEDIT 18 0 "f(x)" "6#-%\"fG6#%\"xG" }{TEXT -1 14 " evaluated at " }{XPPEDIT 18 0 "x = a" "6#/%\"xG%\"aG" }{TEXT -1 59 ". In the Newton-Raphson method of root-finding, we assume " }{XPPEDIT 18 0 "x=a" "6#/%\"xG%\"aG" }{TEXT -1 45 " is close to a root, which we will denote as " }{XPPEDIT 18 0 "x[0]" "6#&%\"xG6#\"\"!" }{TEXT -1 46 ". Then the intersection of the line with the " }{TEXT 259 1 "x" } {TEXT -1 20 " axis, at the point " }{XPPEDIT 18 0 "x[1]" "6#&%\"xG6#\" \"\"" }{TEXT -1 20 ", will be closer to " }{XPPEDIT 18 0 "x[0]" "6#&% \"xG6#\"\"!" }{TEXT -1 6 " than " }{TEXT 263 1 "a" }{TEXT -1 56 ". We then find the intersection of the line tangent to " }{XPPEDIT 18 0 "f (x[1])" "6#-%\"fG6#&%\"xG6#\"\"\"" }{TEXT -1 10 " with the " }{TEXT 260 1 "x" }{TEXT -1 15 " axis, call it " }{XPPEDIT 18 0 "x[2]" "6#&%\" xG6#\"\"#" }{TEXT -1 54 ". We iterate in this fashion until at some p oint the " }{TEXT 261 1 "x" }{TEXT -1 30 " axis intersection point, sa y " }{XPPEDIT 18 0 "x[n]" "6#&%\"xG6#%\"nG" }{TEXT -1 48 ", is satisfa ctorily close to to the actual root " }{XPPEDIT 18 0 "x[0]" "6#&%\"xG6 #\"\"!" }{TEXT -1 3 ". " }}{PARA 0 "" 0 "" {TEXT -1 8 " The " } {TEXT 262 1 "y" }{TEXT -1 46 " intercept of the line described by eq. \+ (1) is" }}{PARA 0 "" 0 "" {TEXT -1 6 "(2) " }{XPPEDIT 18 0 "b = f(a) - diff(f(a),x)*a" "6#/%\"bG,&-%\"fG6#%\"aG\"\"\"*&-%%diffG6$-F'6#F)% \"xGF*F)F*!\"\"" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 8 "and the " }{TEXT 264 1 "x" }{TEXT -1 31 " intercept is found by setting " } {XPPEDIT 18 0 "y = 0" "6#/%\"yG\"\"!" }{TEXT -1 11 ", that is, " } {XPPEDIT 18 0 "x = (-b)/diff(f(a),x)" "6#/%\"xG*&,$%\"bG!\"\"\"\"\"-%% diffG6$-%\"fG6#%\"aGF$F(" }{TEXT -1 14 ". Thus, the (" }{XPPEDIT 18 0 "(k+1)" "6#,&%\"kG\"\"\"\"\"\"F%" }{TEXT -1 4 ")th " }{TEXT 265 1 "x " }{TEXT -1 27 " axis intersection point is" }}{PARA 0 "" 0 "" {TEXT -1 7 "(3) " }{XPPEDIT 18 0 "x[k+1] = ( diff(f(x[k]),x)*x[k] - f(x[k ]) ) / diff(f(x[k]),x)" "6#/&%\"xG6#,&%\"kG\"\"\"\"\"\"F)*&,&*&-%%diff G6$-%\"fG6#&F%6#F(F%F)&F%6#F(F)F)-F26#&F%6#F(!\"\"F)-F/6$-F26#&F%6#F(F %F<" }{TEXT -1 7 " or " }{XPPEDIT 18 0 "x[k+1] = x[k] - f(x[k]) / d iff(f(x[k]),x)" "6#/&%\"xG6#,&%\"kG\"\"\"\"\"\"F),&&F%6#F(F)*&-%\"fG6# &F%6#F(F)-%%diffG6$-F06#&F%6#F(F%!\"\"F;" }{TEXT -1 1 " " }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 38 "2. A Symbolic Newton-Raphson Procedure" } }{EXCHG {PARA 0 "" 0 "" {TEXT -1 123 " We are all used to the Newto n-Raphson method in a purely numerical context. What results if we pe rform the iterations " }{TEXT 270 12 "symbolically" }{TEXT -1 151 "? \+ Here is a procedure that symbolically implements the Newton-Raphson me thod of finding roots as embodied in eq. (3). The arguments are the f unction " }{XPPEDIT 18 0 "f(x)" "6#-%\"fG6#%\"xG" }{TEXT -1 21 ", the \+ starting point " }{TEXT 266 1 "a" }{TEXT -1 38 ", the number of iterat ions to perform " }{TEXT 267 1 "n" }{TEXT -1 10 ", and, if " } {XPPEDIT 18 0 "f(x)" "6#-%\"fG6#%\"xG" }{TEXT -1 65 " is not a procedu re, a fourth argument: the independent variable " }{TEXT 268 2 "x." }} }{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "newt := proc( f::algebraic, a:: algebraic, n::posint )\n local k, d, dk, fk, x0, x;\n #\n # If f(x) is a Maple procedure, we require slightly different\n # treatment fo r the derivatives.\n if type(f,procedure) then\n d := D(f);\n \+ x0 := a;\n for k from 1 to n do\n dk := d(x0);\n fk := f( x0);\n if( not type(fk,numeric) and op([2,1],cost(fk/dk)) > 1000 \+ ) then\n x0 := x0 - normal(fk/dk);\n else\n x0 \+ := x0 - factor(fk/dk);\n fi;\n od;\n #\n # The \"normal\" ca se, where f(x) is just an expression.\n else\n if nargs <> 4 then \n ERROR(`: If f(x) is not a procedure, then the `\n \+ .`4th argument must be the independent variable x.`);\n fi;\n if not type(args[4],name) then\n ERROR(`: Fourth argument x must b e a name.`);\n fi;\n x := args[4];\n d := diff( f, x );\n \+ x0 := a;\n for k from 1 to n do\n dk := subs( x=x0, d );\n \+ fk := subs( x=x0, f );\n if( not type(fk,numeric) and op([2,1 ],cost(fk/dk)) > 1000 ) then\n x0 := x0 - normal(fk/dk);\n \+ else\n x0 := x0 - factor(fk/dk);\n fi;\n od;\n fi ;\nend:" "6#>%%newtGR6%'%\"fG%*algebraicG'%\"aGF)'%\"nG%'posintG7(%\"k G%\"dG%#dkG%#fkG%#x0G%\"xG6\"F6@%-%%typeG6$F(%*procedureGC%>F1-%\"DG6# F(>F4F+?(F0\"\"\"\"\"\"F-%%trueGC%>F2-F16#F4>F3-F(6#F4@%34-F96$F3%(num ericG2\"%+5-%#opG6$7$\"\"#\"\"\"-%%costG6#*&F3FDF2!\"\">F4,&F4FD-%'nor malG6#*&F3FDF2FinFin>F4,&F4FD-%'factorG6#*&F3FDF2FinFinC(@$0%&nargsG\" \"%-%&ERRORG6#(%I:~~If~f(x)~is~not~a~procedure,~then~the~G%Q4th~argume nt~must~be~the~independent~variable~x.G@$4-F96$&%%argsG6#\"\"%%%nameG- F\\p6#%E:~~Fourth~argument~x~must~be~a~name.G>F5&Ffp6#\"\"%>F1-%%diffG 6$F(F5>F4F+?(F0\"\"\"FDF-FEC%>F2-%%subsG6$/F5F4F1>F3-F[r6$/F5F4F(@%34- F96$F3FR2\"%+5-FV6$7$\"\"#\"\"\"-Ffn6#*&F3FDF2Fin>F4,&F4FD-F]o6#*&F3FD F2FinFin>F4,&F4FD-Fco6#*&F3FDF2FinFinF6F6F6" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 11 "3. Examples" }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 25 "3. 1. Roots of a Quadratic" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 " Fir st, let's consider a simple quadratic:" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "f := convert( [seq(a[k]*x^k,k=0..2)], `+` );" "6#>%\"fG -%(convertG6$7#-%$seqG6$*&&%\"aG6#%\"kG\"\"\")%\"xGF0F1/F0;\"\"!\"\"#% \"+G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fG,(&%\"aG6#\"\"!\"\"\"*&& F'6#F*F*%\"xGF*F**&&F'6#\"\"#F*)F.F2\"\"\"F*" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 69 "The second-order Newton-Raphson approximation to one of the roots of " }{TEXT 271 1 "f" }{TEXT -1 22 ", starting at a point \+ " }{XPPEDIT 18 0 "x=b" "6#/%\"xG%\"bG" }{TEXT -1 18 ", is then given b y" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "newt( f, b, 2, x );" "6#-% %newtG6&%\"fG%\"bG\"\"#%\"xG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,(%\"b G\"\"\"*&,(&%\"aG6#\"\"!F%*&&F)6#F%F%F$F%F%*&&F)6#\"\"#F%)F$F2\"\"\"F% F4,&F-F%*&F0F4F$F4F2!\"\"!\"\"*&*&F0F4)F'F2F4F4*&,**$)F-F2F4F%*(F-F4F0 F4F$F4F2*&)F0F2F4F3F4F2*&F0F4F(F%!\"#\"\"\"F5\"\"\"F7F8" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 41 "We could put this in a simpler form, say \+ " }{XPPEDIT 19 1 "factor(%);" "6#-%'factorG6#%\"%G" }{TEXT -1 1 " " }} {PARA 11 "" 1 "" {XPPMATH 20 "6#,$*&,,*&)&%\"aG6#\"\"#\"\"$\"\"\")%\"b G\"\"%F-\"\"\"**&F)6#\"\"!F1&F)6#F1F1F(F1F/F1!\"%*(F3F-)F(F+F-)F/F+F-! \"'*&F3F-)F6F+F-!\"\"*&F(F-)F3F+F-F1F-*&,**$F>F-F?*(F6F-F(F-F/F-!\"#*& F:F-F;F-FF*&F(F-F3F-F+\"\"\",&F6F1*&F(F-F/F-F+\"\"\"!\"\"F?" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 80 "which involves significantly fewer operations, especially for larger iterations " }{TEXT 273 1 "n" } {TEXT -1 117 ". But the former version has the advantage of explicitl y showing the correction terms, making them easy to isolate. " }}}} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 21 "3.2. Roots of a Cubic" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 " Let's consider now a simple cubic, fo r which we know the roots." }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "f := (x-3)*(x+2)*(x-1);" "6#>%\"fG*(,&%\"xG\"\"\"\"\"$!\"\"F(,&F'F(\"\" #F(F(,&F'F(\"\"\"F*F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fG*(,&%\" xG\"\"\"!\"$F(F(,&F'F(\"\"#F(F(,&F'F(!\"\"F(F(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "The second-order Newton-Raphson root approximation, \+ again starting at a point " }{XPPEDIT 18 0 "x=b" "6#/%\"xG%\"bG" } {TEXT -1 4 ", is" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "newt( f, b, 2, x );" "6#-%%newtG6&%\"fG%\"bG\"\"#%\"xG" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#,(%\"bG\"\"\"*&*(,&F$F%!\"$F%F%,&F$F%\"\"#F%F%,&F$F%!\" \"F%F%\"\"\",(*$)F$F+F.\"\"$F$!\"%!\"&F%!\"\"F-*&*.,&F$F+F-F%F%)F,F+F. ,&F$F+F%F%F%)F(F+F.,&F$F%!\"#F%F%)F*F+F.F.*&,0*$)F$\"\"'F.\"#7*$)F$\" \"&F.!#[*$)F$\"\"%F.\"#B*$)F$F2F.\"#cF0\"$u\"F$!$'H!$P\"F%\"\"\"F/\"\" \"F5F=" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "We could also use a Map le procedure for the function " }{XPPEDIT 18 0 "f(x)" "6#-%\"fG6#%\"xG " }{TEXT -1 1 ":" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "f := (x)->( x-3)*(x+2)*(x-1);" "6#>%\"fGR6#%\"xG7\"6$%)operatorG%&arrowG6\"*(,&F' \"\"\"\"\"$!\"\"F/,&F'F/\"\"#F/F/,&F'F/\"\"\"F1F/F,F,F," }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%\"fGR6#%\"xG6\"6$%)operatorG%&arrowGF(*(,&9$\" \"\"!\"$F/F/,&F.F/\"\"#F/F/,&F.F/!\"\"F/F/F(F(F(" }}}{EXCHG {PARA 0 " " 0 "" {XPPEDIT 19 1 "newt( f, b, 2 );" "6#-%%newtG6%%\"fG%\"bG\"\"#" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#,(%\"bG\"\"\"*&*(,&F$F%!\"$F%F%,&F$F %\"\"#F%F%,&F$F%!\"\"F%F%\"\"\",(*$)F$F+F.\"\"$F$!\"%!\"&F%!\"\"F-*&*. ,&F$F+F-F%F%)F,F+F.,&F$F+F%F%F%)F(F+F.,&F$F%!\"#F%F%)F*F+F.F.*&,0*$)F$ \"\"'F.\"#7*$)F$\"\"&F.!#[*$)F$\"\"%F.\"#B*$)F$F2F.\"#cF0\"$u\"F$!$'H! $P\"F%\"\"\"F/\"\"\"F5F=" }}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 35 "4. An Illustrative Convergence Plot" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 " Let's plot our simple test cubic:" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "plot( [f(x),0], x=-3..4, color=[blue,black], labels=[\" x\",\"f(x)\"] );" "6#-%%plotG6&7$-%\"fG6#%\"xG\"\"!/F*;,$\"\"$!\"\"\" \"%/%&colorG7$%%blueG%&blackG/%'labelsG7$Q\"x6\"Q%f(x)F;" }}{PARA 13 " " 1 "" {GLPLOT2D 508 308 308 {PLOTDATA 2 "6+-%'CURVESG6$7U7$$!\"$\"\"! $!#CF*7$$!1nmTg*4P#H!#:$!1B#p/lrp9#!#97$$!1ML$3#*>u%GF0$!1C%REE\"[1>F3 7$$!1+DJl./\"y#F0$!1;%G18Dsq\"F37$$!1m;z43m9FF0$!15#\\X'Q3<:F37$$!1LLe /$f`c#F0$!1'*oe:X\"=7\"F37$$!1LL3K\"o]T#F0$!1>HM)*HyvwF07$$!1m;Hn7\\lA F0$!1RM;Gs'\\c%F07$$!1L$ekO9o7#F0$!1S%p!z*4H.#F07$$!1+]7oCA$)>F0$\"1() o^h\"eT\\#!#;7$$!1L$e9_>Z$=F0$\"1CqZd6=lAF07$$!1+]iDGp'o\"F0$\"1>Y$y>r ]%RF07$$!1mm;u\"HW`\"F0$\"1chu'>G/N&F07$$!1LL3n_J+9F0$\"1&**HIXcRL'F07 $$!1++]sZL\\7F0$\"1ynm42*\\<(F07$$!1++]PVt(4\"F0$\"1rV!HwPev(F07$$!1** ***\\F#R;&*Ffn$\"1s]'y#39&3)F07$$!1KLe9(3(*=)Ffn$\"1j<[MJ8/#)F07$$!1lm mTO:7mFfn$\"1xq%39yD9)F07$$!1jmmmvvv_Ffn$\"1mI.0IOMzF07$$!1***\\7B67s$ Ffn$\"1&Q\">1#G@`(F07$$!1nmm;V$\\\")p[F07$$\"1OL3xplzMFfn$\"1E#*o))H9gSF07$$\"1mm;H([a'\\Ffn$\" 1!GsS`)eYJF07$$\"1omT&yo(3lFfn$\"1UL')G@2u@F07$$\"1-+DJMB_yFfn$\"18M`= :)[K\"F07$$\"1qm;a?@.$*Ffn$\"1h.,?g*eA%Ffn7$$\"1+++&\\@-3\"F0$!1**o]s@ xVZFfn7$$\"1++v$opoA\"F0$!1m0t*Hr!)H\"F07$$\"1,](oMf(o8F0$!11g&\\XFk-# F07$$\"1++DEOIE:F0$!1!*\\q5D/NFF07$$\"1LLLoT'ym\"F0$!1DNtFmCjKF07$$\"1 ,+]i-,>=F0$!1Emga)4Rp$F07$$\"1m;a)3rf&>F0$!1_HRgAI[RF07$$\"1,+]Zaq0@F0 $!1rk=m(H)fSF07$$\"1L$3-\"QfYAF0$!1(R%yzZO))RF07$$\"1+]PWF'QR#F0$!1Uf_ 8AD7PF07$$\"1LL$e/Xy`#F0$!1S_V8&\\^A$F07$$\"1+](=<\"e)o#F0$!1Fv$))R=bY #F07$$\"1ommwzvLGF0$!1$>\\QpfNZ\"F07$$\"1nm\"zAAA)HF0$!1[*=A&*4dv\"Ffn 7$$\"1M$3-7d%HJF0$\"1elm9[099F07$$\"1,++&p]ZE$F0$\"1Yhe:YrcJF07$$\"1ML e*R7)>MF0$\"1))z$z\"*4e]&F07$$\"1nmmO9]eNF0$\"1V6*4F\"pUzF07$$\"1+](o( GP1PF0$\"1T\\L$4#*34\"F37$$\"1,]78Z!z%QF0$\"1'eg+'Q779F37$$\"\"%F*$\"# =F*-%'COLOURG6&%$RGBGF*F*$\"*++++\"!\")-F$6$7S7$F(F*7$F5F*7$F?F*7$FDF* 7$FIF*7$FNF*7$FSF*7$FXF*7$FhnF*7$F]oF*7$FboF*7$FgoF*7$F\\pF*7$FapF*7$F fpF*7$F[qF*7$F`qF*7$FeqF*7$FjqF*7$F_rF*7$FdrF*7$FjrF*7$F_sF*7$FdsF*7$F isF*7$F^tF*7$FctF*7$FhtF*7$F]uF*7$FbuF*7$FguF*7$F\\vF*7$FavF*7$FfvF*7$ F[wF*7$F`wF*7$FewF*7$FjwF*7$F_xF*7$FdxF*7$FixF*7$F^yF*7$FcyF*7$FhyF*7$ F]zF*7$FbzF*7$FgzF*7$F\\[lF*7$Fa[lF*-Ff[l6&Fh[lF*F*F*-%(SCALINGG6#%.UN CONSTRAINEDGFe[l-%+AXESLABELSG6%Q\"x6\"Q%f(x)Fj_l-%%FONTG6%%&TIMESG%&R OMANGFd[l-%*AXESTICKSG6%%(DEFAULTGFd`l-F]`l6$%*HELVETICAG\"#7-%&TITLEG 6$%!G-F]`l6%F_`l%'ITALICGFd[l-%*AXESSTYLEG6#%$BOXG-%%VIEWG6$;F(Fa[lFd` l" 1 2 0 1 0 2 6 1 2 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 140 " Now let's graphically illustrate the convergence to the roots, given several starting points. First defin e a few convenient procedures." }}}{EXCHG {PARA 15 "" 0 "" {TEXT -1 87 "A procedure to create a list of point pairs of the form [iteration number, x intercept]" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "p := p roc(N,b) local n; [ seq( [n,newt(f,b,n)], n=1..N ) ]: end:" "6#>%\"pGR 6$%\"NG%\"bG7#%\"nG6\"F+7#-%$seqG6$7$F*-%%newtG6%%\"fGF(F*/F*;\"\"\"F' F+F+F+" }}}{EXCHG {PARA 0 "" 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 15 "" 0 "" {TEXT -1 92 "A procedure to plot the point pairs as circles and to connect the points with line segments:" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "P := proc( N, b, c )\n global p1, p2;\n p1 := pl ot( p(N,b), style=point, symbol=circle, color=black );\n p2 := plot( \+ p(N,b), style=line, color=c );\n p1,p2;\nend:" "6#>%\"PGR6%%\"NG%\"bG %\"cG7\"6\"F+C%>%#p1G-%%plotG6&-%\"pG6$F'F(/%&styleG%&pointG/%'symbolG %'circleG/%&colorG%&blackG>%#p2G-F06%-F36$F'F(/F6%%lineG/Fplots[display]( \n [ plot([3,-2,1], color=black), P(N,a1,blue), P(N,a2,blue), \n P(N,a3,red), \+ P(N,a4,red), P(N,a5,red) ], \n view=[1..N,-5..7], title=`x i ntercept vs. iteration`, \n titlefont=[HELVETICA,16], axes=b ox ):" "6#>%%doitGR6(%\"NG%#a1G%#a2G%#a3G%#a4G%#a5G7\"6$%)operatorG%&a rrowG6\"-&%&plotsG6#%(displayG6'7(-%%plotG6$7%\"\"$,$\"\"#!\"\"\"\"\"/ %&colorG%&blackG-%\"PG6%F'F(%%blueG-FF6%F'F)FH-FF6%F'F*%$redG-FF6%F'F+ FM-FF6%F'F,FM/%%viewG7$;\"\"\"F';,$\"\"&F@\"\"(/%&titleG%:x~intercept~ vs.~iterationG/%*titlefontG7$%*HELVETICAG\"#;/%%axesG%$boxGF1F1F1" }}} {EXCHG {PARA 0 "" 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Okay. Now create the plot." }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "doit(8,2.2,-1,-0.338,-0.3375,-0.3393);" "6#-%%doitG6(\" \")$\"#A!\"\",$\"\"\"!\"\",$$\"$Q$!\"$F,,$$\"%vL!\"%F,,$$\"%$R$!\"%F, " }}{PARA 13 "" 1 "" {GLPLOT2D 400 300 300 {PLOTDATA 2 "66-%'CURVESG6$ 7S7$$!#5\"\"!$\"\"$F*7$$!1nmm;p0k&*!#:F+7$$!1LL$3s%HaF0F+7$$!1******\\$*4)*\\F0F+7$$!1+++]_&\\c%F0F+7$$!1+++ ]1aZTF0F+7$$!1mm;/#)[oPF0F+7$$!1LLL$=exJ$F0F+7$$!1LLLL2$f$HF0F+7$$!1++ ]PYx\"\\#F0F+7$$!1MLLL7i)4#F0F+7$$!1****\\P'psm\"F0F+7$$!1****\\74_c7F 0F+7$$!1:LL$3x%z#)!#;F+7$$!1ILL3s$QM%FjoF+7$$!1^omm;zr)*!#=F+7$$\"1WLL ezw5VFjoF+7$$\"1.++v$Q#\\\")FjoF+7$$\"1NLLe\"*[H7F0F+7$$\"1++++dxd;F0F +7$$\"1,++D0xw?F0F+7$$\"1,+]i&p@[#F0F+7$$\"1+++vgHKHF0F+7$$\"1lmmmZvOL F0F+7$$\"1,++]2goPF0F+7$$\"1KL$eR<*fTF0F+7$$\"1-++])Hxe%F0F+7$$\"1mm;H !o-*\\F0F+7$$\"1,+]7k.6aF0F+7$$\"1mmm;WTAeF0F+7$$\"1****\\i!*3`iF0F+7$ $\"1NLLL*zym'F0F+7$$\"1OLL3N1#4(F0F+7$$\"1pm;HYt7vF0F+7$$\"1-+++xG**yF 0F+7$$\"1qmmT6KU$)F0F+7$$\"1OLLLbdQ()F0F+7$$\"1++]i`1h\"*F0F+7$$\"1-+] P?Wl&*F0F+7$$\"#5F*F+-%'COLOURG6&%$RGBGF*F*F*-F$6$7S7$F($!\"#F*7$F.Fbu 7$F2Fbu7$F5Fbu7$F8Fbu7$F;Fbu7$F>Fbu7$FAFbu7$FDFbu7$FGFbu7$FJFbu7$FMFbu 7$FPFbu7$FSFbu7$FVFbu7$FYFbu7$FfnFbu7$FinFbu7$F\\oFbu7$F_oFbu7$FboFbu7 $FeoFbu7$FhoFbu7$F\\pFbu7$F_pFbu7$FcpFbu7$FfpFbu7$FipFbu7$F\\qFbu7$F_q Fbu7$FbqFbu7$FeqFbu7$FhqFbu7$F[rFbu7$F^rFbu7$FarFbu7$FdrFbu7$FgrFbu7$F jrFbu7$F]sFbu7$F`sFbu7$FcsFbu7$FfsFbu7$FisFbu7$F\\tFbu7$F_tFbu7$FbtFbu 7$FetFbu7$FhtFbuFjt-F$6$7S7$F($\"\"\"F*7$F.Fhx7$F2Fhx7$F5Fhx7$F8Fhx7$F ;Fhx7$F>Fhx7$FAFhx7$FDFhx7$FGFhx7$FJFhx7$FMFhx7$FPFhx7$FSFhx7$FVFhx7$F YFhx7$FfnFhx7$FinFhx7$F\\oFhx7$F_oFhx7$FboFhx7$FeoFhx7$FhoFhx7$F\\pFhx 7$F_pFhx7$FcpFhx7$FfpFhx7$FipFhx7$F\\qFhx7$F_qFhx7$FbqFhx7$FeqFhx7$Fhq Fhx7$F[rFhx7$F^rFhx7$FarFhx7$FdrFhx7$FgrFhx7$FjrFhx7$F]sFhx7$F`sFhx7$F csFhx7$FfsFhx7$FisFhx7$F\\tFhx7$F_tFhx7$FbtFhx7$FetFhx7$FhtFhxFjt-F$6& 7*7$Fhx$\"1+++++++yF07$$\"\"#F*$\"1+++qt)Qh&F07$F+$\"1+++]!R_C%F07$$\" \"%F*$\"1+++8xleMF07$$\"\"&F*$\"1+++pCn(4$F07$$\"\"'F*$\"1+++2M-1IF07$ $\"\"(F*$\"1+++G_-+IF07$$\"\")F*F+Fjt-%&STYLEG6#%&POINTG-%'SYMBOLG6#%' CIRCLEG-F$6%F\\\\l-F[u6&F]uF*F*$\"*++++\"!\")-F`^l6#%%LINEG-F$6&7*7$Fh x$!\"&F*7$Fa\\l$!1+++++++MF07$F+$!151()pU4*[#F07$Fi\\l$!1zK(y/C74#F07$ F^]l$!1#yD7rIT+#F07$Fc]l$!1y?bp!4++#F07$Fh]l$!1qQ/++++?F07$%%FAILGFj`l FjtF_^lFc^l-F$6%Fc_lFi^lF^_l-F$6&7*7$Fhx$\"1+++4yx2>F07$Fa\\l$!1+++SBU ]NFjo7$F+$\"1+++p[s!)>F07$Fi\\l$!1+++VqAq9F07$F^]l$!1+++/#fWE#F07$Fc]l $!1+++$Q#oI?F07$Fh]l$!1+++8)*[+?F07$F]^l$!1+++G,++?F0FjtF_^lFc^l-F$6%F _al-F[u6&F]uF[_lF*F*F^_l-F$6&7*7$Fhx$\"1+++PDt0>F07$Fa\\l$!1++++O<^LFj o7$F+$\"1+++ps0'*=F07$Fi\\l$!1++++mpgCFjo7$F^]l$\"1+++=SF/;F07$Fc]l$\" 1+++(om!=yFjo7$Fh]l$\"1+++d[Od**Fjo7$F]^l$\"1+++0+(*****FjoFjtF_^lFc^l -F$6%F^clFjblF^_l-F$6&7*7$Fhx$\"1+++W*>J\">F07$Fa\\l$!1+++?CZ!4%Fjo7$F +$\"1+++UdIhAF07$Fi\\l$\"1+++*HClK&F07$F^]l$\"1+++Q[_sSF07$Fc]l$\"1+++ ytapLF07$Fh]l$\"1+++/=#y1$F07$F]^l$\"1+++_Q*H+$F0FjtF_^lFc^l-F$6%F[elF jblF^_l-%(SCALINGG6#%.UNCONSTRAINEDGFi^l-%+AXESLABELSG6%%!GF]gl-%%FONT G6%%&TIMESG%&ROMANG\"#=-%*AXESTICKSG6%%(DEFAULTGFggl-F_gl6$%*HELVETICA G\"#7-%&TITLEG6$%:x~intercept~vs.~iterationG-F_gl6$Fjgl\"#;-%*AXESSTYL EG6#%$BOXG-%%VIEWG6$;FhxF]^l;Fe_lFh]l" 1 2 0 1 0 2 6 1 2 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 387 " The locations of the roots are shown by the black lines. The iteration s tarting points were purposely chosen to be perverse \"guesses\", leadi ng to an artificial need for a higher number of iterations. In practi ce, one would be a little more careful to start close to a root. Noti ce (red curves) that the particular solution settled into can sensitiv ely depend on the initial guess." }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 17 "5. Equation Bloat" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 " By 5th order or so, the " }{TEXT 269 7 "general" }{TEXT -1 185 " solutio ns (i.e., pure symbolics and no numerics) are getting somewhat ugly, e ven for this simple cubic test case. The \"cost\" of the 2nd order so lution shown in the Examples section is" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "cost( newt(f,b,2) );" "6#-%%costG6#-%%newtG6%%\"fG%\"bG \"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,(%*additionsG\"#@%0multiplic ationsG\"#S%*divisionsG\"\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 39 "T he \"cost\" of the 5th order solution is" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "cost( newt(f,b,5) );" "6#-%%costG6#-%%newtG6%%\"fG%\"bG \"\"&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,(%*additionsG\"$q'%0multipli cationsG\"&zp#%*divisionsG\"#:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 184 " Let's make a plot of the log of the total number of additions , multiplications, and divisions as a function of number of iterations . First we write a procedure to take care of it." }}}{EXCHG {PARA 0 " " 0 "" {XPPEDIT 19 1 "costplot := proc( N::posint )\n local p, i, j, \+ c, d, plt1, plt2, ylab, fsiz, \n ymin, ymax;\n p := []:\n for i from 1 to N do\n c := cost( newt(f,b,i) );\n d := 0;\n for j from 1 to nops(c) do\n if nops(op(j,c))=1 then\n d := d + 1;\n else\n d := d + op([j,1],c);\n fi;\n od:\n p := [op(p),[i,log10(d)]];\n od:\n plt1 := plot( p, color=blue ) ;\n plt2 := plot( p, color=red, style=point, symbol=box );\n fsiz := 14;\n ymin := trunc( p[1][2] );\n ymax := trunc( p[nops(p)][2] ) + \+ 1;\n ylab := plots[textplot]( [1.0,ymax-0.5,`log(operations)`], \n \+ font=[HELVETICA,fsiz], align=RIGHT );\n plots [display]( [plt1,plt2,ylab], view=[0.8..N+0.2,ymin..ymax], \n \+ axesfont=[HELVETICA,fsiz-2], axes=box );\nend:" "6#>%)costplo tGR6#'%\"NG%'posintG7-%\"pG%\"iG%\"jG%\"cG%\"dG%%plt1G%%plt2G%%ylabG%% fsizG%%yminG%%ymaxG6\"F6C+>F+7\"?(F,\"\"\"\"\"\"F(%%trueGC&>F.-%%costG 6#-%%newtG6%%\"fG%\"bGF,>F/\"\"!?(F-\"\"\"F<-%%nopsG6#F.F=@%/-FM6#-%#o pG6$F-F.\"\"\">F/,&F/F<\"\"\"F<>F/,&F/F<-FT6$7$F-\"\"\"F.F<>F+7$-FT6#F +7$F,-%&log10G6#F/>F0-%%plotG6$F+/%&colorG%%blueG>F1-Fdo6&F+/Fgo%$redG /%&styleG%&pointG/%'symbolG%$boxG>F3\"#9>F4-%&truncG6#&&F+6#\"\"\"6#\" \"#>F5,&-Fhp6#&&F+6#-FM6#F+6#\"\"#F<\"\"\"F<>F2-&%&plotsG6#%)textplotG 6%7%$\"#5!\"\",&F5F<$\"\"&!\"\"!\"\"%0log(operations)G/%%fontG7$%*HELV ETICAGF3/%&alignG%&RIGHTG-&F_r6#%(displayG6&7%F0F1F2/%%viewG7$;$\"\")! \"\",&F(F<$\"\"#!\"\"F<;F4F5/%)axesfontG7$F`s,&F3F<\"\"#F[s/%%axesGFcp F6F6F6" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 35 "Now create the equation bloat plot." }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "costplot(6);" " 6#-%)costplotG6#\"\"'" }}{PARA 13 "" 1 "" {GLPLOT2D 400 300 300 {PLOTDATA 2 "6,-%'CURVESG6$7(7$$\"\"\"\"\"!$\"1Dw/Y7=z5!#:7$$\"\"#F*$ \"1()Q)R(*zh!=F-7$$\"\"$F*$\"1#3.#)p*>IDF-7$$\"\"%F*$\"1$H=P7OiO$F-7$$ \"\"&F*$\"1[)Hf(\\\">W%F-7$$\"\"'F*$\"1dZ3\"R!=+aF--%'COLOURG6&%$RGBGF *F*$\"*++++\"!\")-F$6&F&-FH6&FJFKF*F*-%&STYLEG6#%&POINTG-%'SYMBOLG6#%$ BOXG-%%TEXTG6'7$$\"#5!\"\"$\"#bFjn%0log(operations)G%+ALIGNRIGHTGFG-%% FONTG6$%*HELVETICAG\"#9-%(SCALINGG6#%.UNCONSTRAINEDGFG-%+AXESLABELSG6% %!GF[p-F`o6%%&TIMESG%&ROMANG\"#=-%*AXESTICKSG6%%(DEFAULTGFdp-F`o6$Fbo \"#7-%&TITLEG6$F[p-F`o6%F^p%'ITALICGF`p-%*AXESSTYLEGFX-%%VIEWG6$;$\"\" )Fjn$\"#iFjn;F(FC" 1 2 0 1 0 2 6 1 2 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 69 "We see that the equation bloat of the purely symbolic solution grows " }{TEXT 272 13 "exponentially " }{TEXT -1 98 " with the number of iterations. In fact, the slope of the curve is approximately (but not quite) " }{XPPEDIT 18 0 "m=1" "6# /%\"mG\"\"\"" }{TEXT -1 85 ", so that with each iteration the operatio n count mushrooms by a factor of almost 10." }}}{EXCHG {PARA 0 "" 0 " " {MPLTEXT 1 0 0 "" }}}}}{MARK "0 0" 0 }{VIEWOPTS 1 0 0 1 1 1803 }