1 / 28

Math.375 Fall 2005 4. Errors in Numerical Computation

Math.375 Fall 2005 4. Errors in Numerical Computation. Vageli Coutsias. Types of errors in numerical computation. Roundoff errors Pi = 3.14159 Pi = 3.1415926535897932384626 Truncation errors Cos x = 1 – x^2/2 Cos x = 1 – x^2/2 + x^4/4! Errors usually accumulate randomly

emanuels
Download Presentation

Math.375 Fall 2005 4. Errors in Numerical Computation

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Math.375 Fall 20054. Errors in Numerical Computation Vageli Coutsias

  2. Types of errors in numerical computation • Roundoff errors Pi = 3.14159 Pi = 3.1415926535897932384626 • Truncation errors Cos x = 1 – x^2/2 Cos x = 1 – x^2/2 + x^4/4! Errors usually accumulate randomly (random walk) But they can also be systematic, and the reasons may be subtle!

  3. x = linspace(1-2*10^-8,1+2*10^-8,401); f = x.^7-7*x.^6+21*x.^5-35*x.^4+35*x.^3-21*x.^2+7*x-1; plot(x,f) CANCELLATION ERRORS

  4. g = -1+x.*(7+x.*(-21+x.*(35+x.*(-35+x.*(21+x.*(-7+x)))))); plot(x,g)

  5. h = (x-1).^7; plot(x,h)

  6. plot(x,h,x,g,x,f)

  7. z(1) = 0; z(2) = 2; for k = 2:29 z(k+1) = 2^(k-1/2)*(1-(1-4^(1-k)*z(k)^2)^(1/2))^(1/2); end semilogy(1:30,abs(z-pi)/pi)

  8. % Plot derivative error vs. h % to exhibit optimal step clear fname = 'sin'; fname1= 'cos'; delta = eps; hmin = 10^6*eps; hmax = 10^10*eps; xmax = 1; x=pi/4; M2=1; % estimate for second derivative hopt = 2*sqrt(delta/M2);% h for minimum error yopt = (M2/2)*hopt+2*delta/hopt; d = (feval(fname,x+hopt)-feval(fname,x))/hopt;

  9. h = linspace(hmin,hmax,1000); y1 = (M2/2)*h; % truncation error y2 = 2*delta./h; % roundoff error y = y1 + y2; % total error derivh =(feval(fname,x+h) - feval(fname,x))./h; abserr = abs(feval(fname1,x)-derivh); loglog(h,y,h,y1,h,y2,h,abserr,hopt,yopt,'bo'); xlabel(sprintf('hopt = %3.2d, d%s(%3.2d)/dx=… %3.2d +/-%3.2d',hopt,fname ,x,d,err)) ylabel('abserr, maxerr') title(sprintf('derivative error as function of h … for %s(x)',fname)) axis([hmin hmax 10^(-12) 10^(-4)])

  10. Vocabulary • Consistency – correctness - … • Convergence (rate) • Efficiency (operation counts) • Numerical Stability (error amplification) • Accuracy (relative vs. absolute error) • Roundoff vs. Truncation

  11. Due to limitations in computer arithmetic, need to practice “defensive coding”. • Mathematicsnumerics + implementation • The sequence of computations carried out for the solution of a mathematical problem can be so complex as to require a mathematical analysis of its own correctness. • How can we determine efficiency? Text

  12. Recurrence relations

  13. Instability is inescapable; But one can learn to ride it!

  14. APPENDIX • NUMERICAL MONSTERS: a collection of curious phenomena exhibited by computations at the limit of machine precision

  15. for k=1:50, term = x*term/k; s = s+term;err(k) = abs(f(k) - s);end relerr = err/exp(x); semilogy(1:nTerms,relerr)% semilogy(1:nTerms,err) ylabel('Relative Error in Partial Sum.') xlabel('Order of Partial Sum.') title(sprintf('x = %5.2f',x)) figure semilogy(1:nTerms,err)%end term = 1; s = 1; f = exp(x)*ones(nTerms,1); for k=1:50, term = x*term/k; s = s+term; err(k) = abs(f(k) - s); end relerr = err/exp(x); semilogy(1:nTerms,relerr)% semilogy(1:nTerms,err) ylabel('Relative Error in Partial Sum.') xlabel('Order of Partial Sum.') title(sprintf('x = %5.2f',x)) figure semilogy(1:nTerms,err)%end The first monster: T(x) = e^x * log(1+e^-1) For x>30 a dramatic deviation from the theoretical value 1 is found. For x>36.7 the plotted value drops (incorrectly) to zero Here: x = linspace(0,40,1000); y = exp(x).*log(1+exp(-x));

  16. Pet Monster I

  17. III.MISCELLANEOUS OPERATIONS: (1) Setting ranges for axes, axis([-40 0 0 0.1]) (2) Superimposing plots: plot(x,y,x,exp(x),'o') (3) using "hold on/hold off" (4) Subplots: subplot(m,n,k)

  18. Summary • Roundoff and other errors • Multiple plots

  19. References • C. Essex, M.Davidson, C. Schulzky, “Numerical Monsters”, preprint. • Higham & Higham, Matlab Guide, SIAM • B5 Trailer; http://www.scifi.com/b5rangers/

More Related