University of Kentucky

MA/CS 321:  Introduction to Numerical Methods

Dr. Qiang Ye

735 Patterson Office Tower,  Phone:  257-4653,  Email:   qye@ms.uky.edu

Interactive Demonstration and Exercises

• This web page contains Octave programs that I have developed for interactive demonstrations and exercises. Click on "Submit to Octave" to run the programs. You can change various parameters in the codes to experiment with the numerical methods. Note that the lines starting with % are comments.
  
  The web-based interface to Octave used here is provided courtesy of Prof. Mai Zhou of Dept. of Statistics.
  
  When a figure does not show on your browser, try by clicking the "refresh" button of the browser.

• Exercise 1: The following plot the graphes of sin(x) and its Taylor polynomials. You can modify the number of terms used in the Taylor polynomial (n) and the range of approximation (a and b) to see how they affect the approximation.
%set the range of plot: a and b; a=-5; b=5; % set the number of terms in Taylor polynomial n n=3 % set ploterror=1 if you only want to plot the error curve ploterror=0; x=(a:0.05:b)'; % set points in x-axis y=sin(x); z=0*x; for i=1:n z=z+(-1)^(i+1)*(x.^(2*i-1))/gamma(2*i); end if (ploterror ==0) data = [x, y, z]; gplot [a:b] [-2:2] data using 1:2, data using 1:3 endif if (ploterror==1); z=y-z; data = [x z]; gplot data endif

• Exercise 2: The following involves subtraction of two close numbers a and b. You can change b by adding more 0s to make it closer to a (e.g. b=100000/500001) and see what happen. Don't forget to change exact_c as well.
format long e a=1/5; b=1000/5001; c=a-b; fprintf('\n a \t= \t%2.20f\n b \t= \t%2.20f\n c \t= \t%2.20f\n\n',a,b,c); % the exact value of c is 1/25005 exact_c=1/25005; fprintf('\n OR IN FLOATING POINT FORM \n c \t\t= \t%1.16e\n exact_c \t= \t%1.16e\n\n',c,exact_c); % relative error rel_error=(c-exact_c)/exact_c; fprintf('\n rel_error \t= \t%1.2e\n',rel_error);

• Exercise 3: The following implements Newton's method for 1/(x^2+1 )-1/2=0. Experiment with
  different initial value x and the number of iterations to observe convergence behavior. The graph of the function is given at bottom. It's also worth trying different equations. Try sin(x) - x =0.
% set initial apprximation here, change this value for experiments x=2; % 20 Newton's iterations fprintf(1,' \nIteration\tApproximation\tError Estimate\n\n'); for i=1:20, fx=1/(x^2+1)-1/2; fprime=-2*x/(x^2+1)^2; xnew=x-fx/fprime; error=abs(x-xnew); fprintf(1,' %d\t\t%2.16e\t\t%e\n',i,xnew,error); x=xnew; end %plot the function x=(-20:0.1:20)'; y=1./(x.^2+1)-1/2; data = [x y 0*y]; gplot data using 1:2, data using 1:3

• Exercise 4: The following computes the interpolating polynomial for a given input of data points. Experiment different coordinates (x and y) and try more points.
% set the data points for interpolation % x-coordinates in x and y-coordinates in y x=[0; 1; 2; 3; 4]; y=[1; 2; 0; -4; -3]; u=(min(x)-1:0.05:max(x)+1)'; % set points in x-axis % compute the interpolating polynomial n = length(x); v = zeros(size(u)); for k = 1:n w = ones(size(u)); for j = [1:k-1 k+1:n] w = (u-x(j))./(x(k)-x(j)).*w; end v = v + w*y(k); end plot(u,v)

Next interpolation exercise uses points on the graph of sin(x) so that the interpolation approximate sin(x). Experiment different number of interpolation points. Also compare it with Taylor approximation.
% set the data points for interpolation % x-coordinates in x and y-coordinates in y x=0:2:6; x=x'; y=sin(x); u=(min(x)-1:0.05:max(x)+0.1)'; % set points in x-axis % compute the interpolating polynomial n = length(x); v = zeros(size(u)); for k = 1:n w = ones(size(u)); for j = [1:k-1 k+1:n] w = (u-x(j))./(x(k)-x(j)).*w; end v = v + w*y(k); end w=sin(u); data = [u, v, w]; gplot [0:6] [-2:2] data using 1:2, data using 1:3

• Exercise 5: The following implements Richardson extrapolation for approximation f'(x) by
  Method 1. (f(x+h)-f(x))/h
  Method 2. (f(x+h)-f(x-h))/(2h)
  We test it for f(x)=log(1+x) at x=0. Modify the coefficient c used for Richardson extrapolation and observe and explain what value of c makes it work.


Method 1:
format short e % set the coefficient used in Richardson extrapolation c c=2; % initial approximation with h=1 h=1; fprime=(log(1+h)-log(1))/h; error=abs(fprime-1); fprintf(1,' \nh\tError\t\tError of Richardson\n\n') for i=1:10, h=h/2; fprimeold=fprime; fprime=(log(1+h)-log(1))/h; error=abs(fprime-1); richardson=fprime+(fprime-fprimeold)/(c-1); error_richardson=abs(richardson-1); fprintf(1,' %e\t\t%e\t%e\n',h,error,error_richardson); end

Method 2:
format short e % set the coefficient used in Richardson extrapolation c c=4; % initial approximation with h=1 h=1; fprime=(log(1+h)-log(1-h))/(2*h); error=abs(fprime-1); fprintf(1,' \n h \tError \t\tError of Richardson\n\n'); for i=1:10, h=h/2; fprimeold=fprime; fprime=(log(1+h)-log(1-h))/(2*h); error=abs(fprime-1); richardson=fprime+(fprime-fprimeold)/(c-1); error_richardson=abs(richardson-1); fprintf(1,' %e\t\t%e\t%e\n',h,error,error_richardson); end