Newton’s method controlPatrick Dale McCray, Ph.D.29 March, 2011IntroductionThe typical way of controling an iterative scheme is to compare successive values of the quantity being computed with thepreceding value. In order to achieve high accuracy in computed quantities, this notebook introduces a Module based on theidea of having Mathematica perform the iteration with exact values and inspecting the resulting exact values numerically with theN function.ü Exact computationDecimal numbers with a few digits are stored with limited machine precision. Here we compute a value z using (approximate)machine numbers.In[130]:=z 12.4567 37.8; z, Nz,Nz, 100Out[130]=0.329542, 0.329542, 0.329542On the other hand, numbers calculated with whole numbers are stored exactly, such as w.In[131]:=w 124567 10 000378 10; w, Nw,Nw, 100Out[131]=124 567378 000, 0.329542,0.3295423280423280423280423280423280423280423280423280423280423280423280423280423280423280423280423280How do these two (same (??)) values compare?In[132]:=z wOut[132]=5.55112 1017The correct, exact value of the computation, namely 124567/378000, is rational, hence has an eventually periodic decimalexpansion that never terminates, unlike the machine version, 0.329542, which is finite (and incorrect!).Sometimes extra precision may be required. The default value of the maximum extra precision is 50.In[133]:=$MaxExtraPrecisionOut[133]=50We will have to increase this in the following discussion.Code for Newton’s methodSuppose we wish to solve f(x) = 0 for x. The basic Newton step may be computed as follows.In[134]:=NewtonStepx_ : x fxDft,t.t xFor example, evaluate the Newton step at x = 3 would yieldIn[135]:=NewtonStep3Out[135]=3 f3f3Using the function NewtonStep we may write a module for iterating until a desired level of accuracy is obtained. Thefollowing Module is intended to iterate until the iterates agree to within 40 decimal places and takes one argument, the initialvalue for the interation, namely x0.In[136]:=NewtonIterationxinit_ :Module k, test, x0, change local variables ,k 0; test True; x0 xinit;Whiletest is true, repeat the following ,xk NewtonStepx0;k k 1;change Absxk x0;test change 1 10^100;Print"x"k," ", Nxk, 45,", ",Nchange, 9;x0 xk ; Print"Stop"ü ExampleLet us solve the following equation:In[137]:=x 3 SinxOut[137]=x3 Sinx2 Newtons method control.nbIn[138]:=Plotx 3, Sinx, x, 0, PiOut[138]=0.5 1.0 1.5 2.0 2.5 3.00.20.40.60.81.0First we code the Newton step.In[139]:=gx_ : 3 Sinx xIn[140]:=NewtonStepx_ : x gxg'xIn[141]:=NewtonStepxOut[141]=x x 3 Sinx 1 3 CosxNow we invoke our module with an initial value of 2, first checking the default value of the maximum extra precision. In[142]:=$MaxExtraPrecisionOut[142]=50In[143]:=NewtonIteration2Newtons method control.nb 3x1 2.32373206111338325354188208345366581392706335, 0.323732061x2 2.27959482163992639068772151841669254197757185, 0.0441372395x3 2.27886286684763997073119345250115932298600016, 0.000731954792x4 2.27886266007584482041679741117310186574317502, 2.06771795 107x5 2.27886266007582831269995110466710382941208700, 1.65077168 1014x6 2.27886266007582831269995110456188862881827474, 1.05215201 1028N::meprec : Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -2 +-2 + 3 Sin2-1 + 3 Cos2+á1à-1á1à 3 á1àá1à+á3à+-2 +-2 + 3á1à-1 + 3á1à+á3à+á1àá1à+ 3Siná1à-1 + 3 Cos2 +á4à+Timesá3à+ 3Siná1à. àN::meprec : Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating á1à1 - 3Cosá1à. àx7 2.27886266007582831269995110456188862881827474, 4.3 1057N::meprec : Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -2 +-2 + 3 Sin2-1 + 3 Cos2+á1à-1á1à 3 á1àá1à+á4à+-2 +-2 + 3á1à-1 + 3á1à+á4à+á1àá1à+ 3Siná1à-1 + 3 Cos2 +á5à+Timesá3à+ 3Siná1à. àGeneral::stop : Further output of N::meprec will be suppressed during this calculation. àx8 2.27886266007582831269995110456188862881827474, 0. 1059StopSince the internal precision limit has been exceeded, we increase it and try again.In[144]:=$MaxExtraPrecision 200Out[144]=200In[145]:=NewtonIteration2x1 2.32373206111338325354188208345366581392706335, 0.323732061x2 2.27959482163992639068772151841669254197757185, 0.0441372395x3 2.27886286684763997073119345250115932298600016, 0.000731954792x4 2.27886266007584482041679741117310186574317502, 2.06771795 107x5 2.27886266007582831269995110466710382941208700, 1.65077168 1014x6 2.27886266007582831269995110456188862881827474, 1.05215201 1028x7 2.27886266007582831269995110456188862881827474, 4.27426496 1057x8 2.27886266007582831269995110456188862881827474, 7.05386830 10114StopNow we reset the internal precision limit back to its original default value.4 Newtons method control.nbIn[146]:=$MaxExtraPrecision 50Out[146]=50ü Investigating the final iterateSolution achieved for iterates agreeing to withing 100 places:In[147]:=Nxk, 100Out[147]=2.278862660075828312699951104561888628818274740739776516525585529248344464701839186256781340580146512Mathematica high precision calculation:In[148]:=Xk x . FindRootgx 0, x, 2.2, WorkingPrecision 250Out[148]=2.278862660075828312699951104561888628818274740739776516525585529248344464701839186256781340580146512389613012252192308094868943884998837137383845350699117003501335444216009373218571499740602576404645321970635956073298400582545637493233865175964304189The difference:In[149]:=xk XkOut[149]=1.9211406436879307525509 10227As we may see, since the computations are done EXACTLY, the corresponding expressions for x1, x2, x3, get rather LARGE:In[150]:=x1 NewtonStep2Out[150]=2 2 3 Sin2 1 3 Cos2In[151]:=x2 NewtonStepx1Out[151]=2 2 3 Sin2 1 3
View Full Document
Unlocking...