echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 4.2 (Polynomial Calculus).
% Section 4.2, Introduction to Interpolation, Page 212
echo on; clc; format short; hold off; clear

% Proceed with an investigation of synthetic division for

% evaluating Pn(x), P'n(x) and ºPn(x)dx where

% Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)

% It is assumed that the constant K of integration is K = 0.

% Pn(x), P'n(x) and ºPn(x)dx are evaluated by invoking the Matlab

% subroutine polyval(C,X) polyval(Cd,X) polyval(Ci,X), respectively.

pause % Press any key to continue.

clc;
% Several Taylor polynomial coefficient lists are

% stored in M-files named; zcos zsin ztan zexp

% zacos zasin zatan zcosh zsinh zsqrt zlog

% zsqrt4 zinv zemx2d2 zcosde zsinde zlogq

% The user must determine the appropriate formula

% for f'(x) in order to graph f'(x) and P'n(x).

% The user must also determine the formula for ºf(x)dx

% in order to graph ºf(x)dx and ºPn(x)dx.

pause % Press any key continue.

clc;
% The algorithms are illustrated with the coefficients of the

% Taylor polynomial of degree n = 11 for f(x) = ln(1+x).

% Issue the command zlog to load the coefficients

% into the array C. The function name is loaded

% into the variable fun, the maximum degree is loaded into N.

% For graphing over [a,b] the endpoints are in a and b.

% Load the Taylor coefficients.

[fun,dfun,ifun,x0,m,C,Ax] = zlog;

pause % Press any key continue.

clc;
a = Ax(1,1); % You can change the endpoint a.
b = Ax(1,2); % You can change the endpoint b.
c = Ax(1,3);
d = Ax(1,4);
n = 11; % You can change the degree n.
C = flipud(C);
C = C(1:n+1);
C = flipud(C);

pause % Press any key to graph f(x) and P(x).

clc; clg;
h = (b-a)/200;
X = a:h:b;
x = X;
Y = eval(fun);
Z = polyval(C,X);
axis([a b c d]);...
plot(X,Y,'-g',X,Z,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = [fun,' and P'];...
Mx2 = [Mx1,num2str(n),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'The function is f(x) = ';
Mx2 = 'The interval is ';
Mx3 = 'Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
Mx4 = 'The degree is n = ';
Mx5 = ', and the coefficients list c is:';
clc,echo off,diary output,...
disp(''),disp([Mx1,fun]),...
disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
disp(Mx3),disp([Mx4,num2str(n),Mx5]),...
for i=1:5:n+1, disp(C([i:min(i+4,n+1)])'); end,...
diary off, echo on

pause % Press any key to graph f`(x) and P`n(x).

clc; clg;
a = Ax(2,1);
b = Ax(2,2);
c = Ax(2,3);
d = Ax(2,4);
x = X;
Y = eval(dfun);
Dv = zeros(n,1);
for j=1:n,
Dv(j) = n-j+1; % Form vector [n,n-1,...,1]
end
Cd = C(1:n).*Dv; % Coefficients for Pn`(x).
Z = polyval(Cd,X);
axis([a b c d]);...
plot(X,Y,'-g',X,Z,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = [dfun,' and P`'];...
Mx2 = [Mx1,num2str(n),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'The function is f`(x) = ';
Mx2 = 'The interval is ';
Mx3 = 'P`n(x) = d(1)x^n + d(2)x^(n-1) + ... + d(n)x + d(n+1)';
Mx4 = 'The degree is n = ';
Mx5 = ', and the coefficients list Cd is:';
clc,echo off,diary output,...
disp(''),disp([Mx1,dfun]),...
disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
disp(Mx3),disp([Mx4,num2str(n-1),Mx5]),...
for i=1:5:n, disp(Cd([i:min(i+4,n)])'); end,
diary off, echo on

pause % Press any key to graph ºf(x)dx and ºPn(x)dx.

clc; clg;
a = Ax(3,1);
b = Ax(3,2);
c = Ax(3,3);
d = Ax(3,4);
x = X;
Y = eval(ifun);
Iv = zeros(n+1,1);
for j=1:n+1,
Iv(j) = 1/(n-j+2); % Form vector [1/(n+1),...,1/2,1]
end
Ci = C(1:n+1).*Iv; % Coefficients for ºPn(x)dx.
Ci = [Ci;0];
Z = polyval(Ci,X);
axis(Ax(3,:));...
plot(X,Y,'-g',X,Z,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = [ifun,' and ºP'];...
Mx2 = [Mx1,num2str(n),'(x)dx'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'The function is ºf(x)dx = ';
Mx2 = 'The interval is ';
Mx3 = 'ºPn(x)dx = a(1)x^n + a(2)x^(n-1) + ... + a(n)x + 0';
Mx4 = 'The degree is n = ';
Mx5 = ', and the coefficients list Ci is:';
clc,echo off,diary output,...
disp(''),disp([Mx1,ifun]),...
disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
disp(Mx3),disp([Mx4,num2str(n+1),Mx5]),...
for i=1:5:n+2, disp(Ci([i:min(i+4,n+2)])'); end,
diary off, echo on