C C In this example we experimentally determine the order of convergence C of Newton's method applied to the root problem f(x)=2-exp(x)=0, whose C root is obviously alpha=log(2)=ln(2) C >> help newton NEWTON: [root,vals,iter,ierr]=newton(f,fp,x0,tol,itermax) Implements Newton's Method to find a root of f(x). Input: f - a function of x as a string, e.g., f='x^2-3' fp - the derivative of f as a string, e.g., fp='2*x' x0 - initial guess tol - stopping critera; when abs(f(x_cur)) < tol the routine terminates itermax - maximum number of iterations Output: root - approximate root of f vals - a history of the iterations iter - number of iterations ierr - error flag: 0 - no error 1 - itermax reached >> [root,vals]=newton('2-exp(x)','-exp(x)',5,1e-20,100); x_n f(x_n) ------------------------------------------------ 5.000000000000000e+00 -1.464131591025766e+02 4.013475893998171e+00 -5.333888876400096e+01 3.049616843783578e+00 -1.910725549682328e+01 2.144370991252017e+00 -6.536669908113822e+00 1.378654385307006e+00 -1.969556538483261e+00 8.824890027537665e-01 -4.169079165495138e-01 7.099926020035425e-01 -3.397621124228500e-02 6.932882713164430e-01 -2.822014205330348e-04 6.931471905127781e-01 -1.990566556031581e-08 6.931471805599454e-01 0 >> alpha=log(2) alpha = 6.931471805599453e-01 C C We determine the error in each iterate C >> err=abs(vals(:,1)-alpha) err = 4.306852819440055e+00 3.320328713438226e+00 2.356469663223633e+00 1.451223810692071e+00 6.855072047470608e-01 1.893418221938212e-01 1.684542144359724e-02 1.410907564977082e-04 9.952832780157905e-09 1.110223024625157e-16 C C We need to know the number of iterates, and use the MATLAB command C command 'length'. C >> n=length(err) n = 10 C C Now apply the order formula discussed in class. Note the use C of MATLAB vector operations using a '.' before the operator. C >> log(err(3:n)./err(2:n-1))./log(err(2:n-1)./err(1:n-2)) ans = 1.318114822613668e+00 1.413702317303563e+00 1.547172873725818e+00 1.715465438359840e+00 1.880511101149659e+00 1.976639845998064e+00 1.998837360536195e+00 1.915557544776791e+00 C C Note that the order appears to be converging to 2 except for the C last iterate. So we take the order=2, just as the theory predicts C since alpha=log(2) is a SIMPLE root! C >> order=2; C C To determine the constant C apply C is approximately C (e_{n+1}/(e_{n})^order in the limit. C >> err(2:n)./(err(1:n-1).^order) ans = 1.790032927495308e-01 2.137468340883087e-01 2.613428899066555e-01 3.254941839662007e-01 4.029234140199905e-01 4.698819647181383e-01 4.972042137121797e-01 4.999764899906062e-01 1.120770819036981e+00 C C Again we ignore the last iterate. C=0.5, which is approximately C abs[0.5*f``(alpha)/f`(alpha)]=0.5, just as the theory predicts!