Hangjie Ji Homework 1

Suppose you are to make a table of values of log(x), 1/2<=x<=3/2, on an equally space mesh of width 1/(N+1), for some positive integer N. Assume that polynomial interpolation will be applied to the table, and that the maximum error in interpolation is to be at most 10^(-6). What value of N should you choose? To answer this question, write a computer program to

Solution: Given N interpolation points (endpoints 1/2 and 3/2 not included), we have N+2 pairs of (x_j,f(x_j)), where f(x)=log(x). To apply the Newton interpolation in this case, first we need to generate the coefficients via divided differences. Let n = N+1 and execute Algorithm 10.2.1, we can find the divided differences \gamma_j=f[x_0,x_1,\cdots,x_j]. Then apply Algorithm 10.2.2, we get the interpolation polynomial p(x).

Set N = 10, Figure 1 shows the polynomial interpolation for the function f(x)=log(x).

Randomly choose a point x0 = 1, the error in the interpolation at x0 is 3.7055e-07.

With N = 10, the error log(x)-p(x) as a function of x is plotted in Figure 2. It is obvious that the error becomes relatively large near the endpoints x = 1/2 or x = 3/2. And the error is quite small when x is away from the endpoints.

Figure 3 shows a semilog plot of the maximum error as a function of N, where the maximum error is evaluated at points x_j=1/2+j/100(N+1) for 0 \le j \le 100(N+1). According to the figure, the maximum error in interpolation is to be at most 10^{-6} if 10<=N<=47.

% Interpolating polynomial p(x) for log(x).
clear all;
close all;

N = 10; %The number of interpolation points
x0 = 1; % arbitrary point for error evaluation
x = linspace(.5,1.5,N+2);
xx = linspace(0.5,1.5,100*(N+1)+1);
n = size(x,2);
f = log(x);
r = f;

for j = 2 : n
    for i = n:-1:j
        r(i)=(r(i)-r(i-1))/(x(i)-x(i-j+1));
    end
end

for i = 1:100*(N+1)+1
    p = r(n);
    for j = n-1:-1:1
        p = r(j)+p*(xx(i)-x(j));
    end
    P(i) = p;
end


figure(1),plot(x,f,'r-o',xx,P,'b-*');
xlabel('x','FontSize',14),ylabel('f(x) & p(x)','FontSize',14);

err0 = log(x0)-P(50*(N+1)+1);
disp('x0'); disp(err0);  % interpolating error for an arbitrary x0.

err = log(xx)-P;
disp(max(err));
figure(2),plot(xx,err,'r-o');
xlabel('x','FontSize',14),ylabel('error','FontSize',14);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nlength = 60;
merr = zeros(Nlength,1);


for N = 1 : Nlength
    x = linspace(.5,1.5,N+2);
    n = size(x,2);
    f = log(x);
    r = f;
    xx = linspace(0.5,1.5,100*(N+1)+1);
    P = zeros(1,100*(N+1)+1);
    for j = 2 : n
        for i = n:-1:j
            r(i)=(r(i)-r(i-1))/(x(i)-x(i-j+1));
        end
    end
    for i = 1:100*(N+1)+1
        p = r(n);
        for j = n-1:-1:1
            p = r(j)+p*(xx(i)-x(j));
        end
        P(i) = p;
    end
    err = log(xx)-P;
    merr(N) = max(err);
end

figure(4), semilogy(1:Nlength,merr,'r-o');
xlabel('N','FontSize',14),ylabel('max(err)','FontSize',14);

Math 563 Scientific Computing II

Homework 1 Due: Jan.15

Hangjie Ji

Figure 1

Figure 2

Figure 3