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.p (Pade rational Approximation).
% Section 4.6, Pade Approximations, Page 249
echo on; clc; format short; hold off; clear

% Investigation of a Pade rational approximations.
% Pn(x)
% f(x) Å Rn,m(x) = -------
% Qm(x)
% where the degrees of Pn and Qm are n and m, respectively.

% The algorithm requires the Maclaurin coefficients of f(x).
% Coefficient lists for several functions have been
% stored in M-files named; zcos zsin ztan zexp
% zacos zasin zatan zcosh zsinh zsqrt zlog
% zsqrt4 zinv zemx2d2 zcosde zsinde zlogq

% Remark. pade.m is used for Algorithm 4.p

pause % Press any key continue.

clc;
% Pade rational approximations for f(x) = cos(x).

% Issue the command zcos to load the coefficients

% into the array C. The function name is loaded

% into the variable fun, the degree is loaded into m.

% The endpoints of [a,b] are loaded into a and b.

% Load the Taylor coefficients.

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

pause % Press any key to continue.

clc;

% Place the degree of Pn(x) in n

% Place the degree of Qm(x) in m

% Remark. Some combinations of n and m are incompatible.
% This might result in a singular matrix
% or division by zero where Qm(x) = 0.

n = 2;
m = 2;
[P,Q] = pade(C,n,m);

pause % Press any key to graph f(x) and Rn,m(x).

clc; clg;
a = -3; % You can change the left endpoint a.
b = 3; % You can change the right endpoint b.
c = -1.5;
d = 1;
X = a:.05:b;
x = X;
Y = eval(fun);
axis([a b c d]);...
R = polyval(P,X)./polyval(Q,X);...
plot(X,Y,'-g',X,R,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['Comparison of ',fun,' and R'];...
Mx2 = [Mx1,num2str(n),',',num2str(m),'(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) = p(1)x^n + p(2)x^(n-1) + ... + p(n)x + p(n+1)';
Mx4 = 'The degree is n = ';
Mx5 = ', and the coefficients list P is:';
Mx6 = 'Qm(x) = q(1)x^m + q(2)x^(m-1) + ... + q(m)x + q(m+1)';
Mx7 = 'The degree is m = ';
Mx8 = ', and the coefficients list Q is:';
clc,format short e,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(P([i:min(i+4,n+1)])); end,...
disp(Mx6),disp([Mx7,num2str(m),Mx8]),...
for i=1:5:m+1, disp(Q([i:min(i+4,m+1)])); end,...
diary off, echo on

pause % Press any key to graph f(x) - Rn,m(x).

clc; clg;
a = -1; % You can change the left endpoint a.
b = 1; % You can change the right endpoint b.
h = (b-a)/200;
X = a:h:b;
x = X;
Y = eval(fun);
R = polyval(P,X)./polyval(Q,X);
c = min(Y-R);
d = max(Y-R);
axis([a b c d]);...
plot(X,Y-R,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['The error; ',fun,' - R'];...
Mx2 = [Mx1,num2str(n),',',num2str(m),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key for a list of numerical computations.

clc; format long;
a = -1.2; % You can change the left endpoint a.
b = 1.2; % You can change the right endpoint b.
X = a:0.1:b;
x = X;
Y = eval(fun);
R = polyval(P,X)./polyval(Q,X);
points = [X;Y;R;Y-R];
Mx1=['Pade rational approximation of f(x) = ',fun];
Mx2=' x(k) f(x(k)) Rn,m(x(k)) error';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
diary off,echo on

pause % Press any key to construct a Pade rational approximation.

clc;

% Place the degree of Pn(x) in n

% Place the degree of Qm(x) in m

% Remark. Some combinations of n and m are incompatible.
% This might result in a singular matrix
% or division by zero where Qm(x) = 0.

n = 4;
m = 4;
[P,Q] = pade(C,n,m);

pause % Press any key to graph f(x) and Rn,m(x).

clc; clg;
a = -5; % You can change the left endpoint a.
b = 5; % You can change the right endpoint b.
c = -1;
d = 1;
X = a:.05:b;
x = X;
Y = eval(fun);
axis([a b c d]);...
R = polyval(P,X)./polyval(Q,X);...
plot(X,Y,'-g',X,R,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['Comparison of ',fun,' and R'];...
Mx2 = [Mx1,num2str(n),',',num2str(m),'(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) = p(1)x^n + p(2)x^(n-1) + ... + p(n)x + p(n+1)';
Mx4 = 'The degree is n = ';
Mx5 = ', and the coefficients list P is:';
Mx6 = 'Qm(x) = q(1)x^m + q(2)x^(m-1) + ... + q(m)x + q(m+1)';
Mx7 = 'The degree is m = ';
Mx8 = ', and the coefficients list Q is:';
clc,format short e,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(P([i:min(i+4,n+1)])); end,...
disp(Mx6),disp([Mx7,num2str(m),Mx8]),...
for i=1:5:m+1, disp(Q([i:min(i+4,m+1)])); end,...
diary off, echo on

pause % Press any key to graph f(x) - Rn,m(x).

clc; clg;
a = -1; % You can change the left endpoint a.
b = 1; % You can change the right endpoint b.
h = (b-a)/200;
X = a:h:b;
x = X;
Y = eval(fun);
R = polyval(P,X)./polyval(Q,X);
c = min(Y-R);
d = max(Y-R);
axis([a b c d]);...
plot(X,Y-R,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['The error; ',fun,' - R'];...
Mx2 = [Mx1,num2str(n),',',num2str(m),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key for a list of numerical computations..

clc; format long;
a = -1.2; % You can change the left endpoint a.
b = 1.2; % You can change the right endpoint b.
X = a:0.1:b;
x = X;
Y = eval(fun);
R = polyval(P,X)./polyval(Q,X);
points = [X;Y;R;Y-R];
Mx1=['Pade rational approximation of f(x) = ',fun];
Mx2=' x(k) f(x(k)) Rn,m(x(k)) error';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
diary off,echo on