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 5.5 (Trigonometric Polynomials).
% Section 5.4, Fourier Series and Trigonometric Polynomials, Page 311
echo on; clc; format short; hold off; clear

% This program finds the trigonometric polynomial

% for periodic data in the interval [-pi,pi].

% Remark. tpcoeff.m and tp.m are used for Algorithm 5.5

pause % Press any key to continue.

clc;
% The computer will place the abscissas in X
% The computer will place the ordinates in Y
% Place the f(x) in the 'string function' fun
% Place the number of points in n.
% Place the order of the trig. poly. in m

fun = 'x/2';
n = 12;
m = 5;

h = 2*pi/n;
X = -pi:h:pi;
x = X;
Y = eval(fun);

[A,B] = tpcoeff(X,Y,m);

pause % Press any key to find the trigonometric polynomial.

clc;
hs = 2*pi/150;
Xs = -pi:hs:pi;
Ys = tp(A,B,Xs,m);
a = -3.5;
b = 3.5;
c = -2;
d = 2;
axis([a b c d]);...
plot(X,Y,'or',Xs,Ys,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = 'Trigonometric polynomial: y = P';...
Mx2 = [Mx1,num2str(m),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;
format short e;
Mx1 = 'The trigonometric polynomial has been determined.';
Mx2 = 'The degree is m = ';
Mx3 = [Mx2,num2str(m)];
Mx4 = 'The coefficients a(n) for cos(nx) are:';
Mx5 = 'The coefficients b(n) for sin(nx) are:';
clc,echo off,diary output,...
disp(Mx1),disp(''),disp(Mx3),disp(''),...
disp(Mx4),disp(''),for i=1:5:m+1,disp(A([i:min(i+4,m+1)]));end,...
disp(Mx5),disp(''),for i=1:5:m+1,disp(B([i:min(i+4,m+1)]));end,...
diary off,echo on

pause % Press any key to continue.

clc;
format short;
points = [X;Y;tp(A,B,X,m);Y-tp(A,B,X,m)]';
Mx6 = ' x(k) y(k) P(x(k)) error';
clc,echo off,diary output,...
disp(''),disp(Mx6),disp(points),diary off,echo on

pause % Press any key to find another trigonometric polynomial.

clc;
% Place the order of the trig. poly. in m

m = 2;

[A,B] = tpcoeff(X,Y,m);

% Remark. Now the degree is only m = 2.

pause % Press any key to continue.

hs = 2*pi/150;
Xs = -pi:hs:pi;
Ys = tp(A,B,Xs,m);
a = -3.5;
b = 3.5;
c = -2;
d = 2;
axis([a b c d]);...
plot(X,Y,'or',Xs,Ys,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = 'Trigonometric polynomial: y = P';...
Mx2 = [Mx1,num2str(m),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;
format short e;
Mx1 = 'The trigonometric polynomial has been determined.';
Mx2 = 'The degree is m = ';
Mx3 = [Mx2,num2str(m)];
Mx4 = 'The coefficients a(n) for cos(nx) are:';
Mx5 = 'The coefficients b(n) for sin(nx) are:';
clc,echo off,diary output,...
disp(Mx1),disp(''),disp(Mx3),disp(''),...
disp(Mx4),disp(''),for i=1:5:m+1,disp(A([i:min(i+4,m+1)]));end,...
disp(Mx5),disp(''),for i=1:5:m+1,disp(B([i:min(i+4,m+1)]));end,...
diary off,echo on

pause % Press any key to continue.

format short;
points = [X;X;tp(A,B,X,m);Y-tp(A,B,X,m)]';
Mx6 = ' x(k) y(k) P(x(k)) error';
clc,echo off,diary output,...
disp(''),disp(Mx6),disp(points),diary off,echo on