%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sheet 2 Q3(a): Truncated cosine and sine series for
% f(x) = exp(x) for 0 <= x <= 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Enter "FS_PDE_Sheet2_Q3;" in command line to run

% truncate after 2*Nt+1 terms 
Nt = 50;

% plot domain
a = -2;
b = 2;

% symbolic variables
syms x n

% where n is a natural number
assume(n,'integer');
assumeAlso(n > 0);

% f_e and f_o as symbolic functions for |x| <= 1
syms fe(x) fo(x)
fe(x) = exp(abs(x));
fo(x) = sign(x)*exp(abs(x));

% evaluate Fourier coefficients
disp('Fourier coefficients');
fe_a0 = simplify( int(fe,-1,1) )
fe_an = simplify( int(fe*cos(n*pi*x),-1,1) )
fe_bn = simplify( int(fe*sin(n*pi*x),-1,1) )
%
fo_a0 = simplify( int(fo,-1,1) )
fo_an = simplify( int(fo*cos(n*pi*x),-1,1) )
fo_bn = simplify( int(fo*sin(n*pi*x),-1,1) )

% construct fourier series as symbolic expressions

Sfe = fe_a0/2;
Sfo = fo_a0/2;

for k=1:Nt
    Sfe = Sfe + subs(fe_an,n,k)*cos(k*pi*x) + subs(fe_bn,n,k)*sin(k*pi*x);
    Sfo = Sfo + subs(fo_an,n,k)*cos(k*pi*x) + subs(fo_bn,n,k)*sin(k*pi*x);
end

% turn it into a matlab function for plotting
mSfe = matlabFunction(Sfe);
mSfo = matlabFunction(Sfo);

% sample points in domain
xx = linspace(a,b,1000);

% sample values of Fourier series
valSfe = mSfe(xx);
valSfo = mSfo(xx);

% sample values of f
% mod operation to sample values outside interval -1 < x <= 1
valfe = fe( mod(xx-1,2)-1 );
valfo = fo( mod(xx-1,2)-1 );

% plot
figure(1);
set(gcf,'Position',[0 0 900 600],'Color','w');
%
subplot(2,1,1)
plot(xx,valSfe,'r-','LineWidth',1.25);
hold on;
grid on;
axis([a b 0 4]);
plot(xx,valfe,'k');
legend([ 'Cosine series, ' num2str(Nt+1) ' terms'],'f_e(x)');
xlabel('x')
ylabel('y')
title('f(x) = exp(x) for 0 <= x <= 1');
%
subplot(2,1,2)
plot(xx,valSfo,'r-','LineWidth',1.25);
hold on;
grid on;
axis([a b -4 4]);
plot(xx,valfo,'k');
legend([ 'Sine series, ' num2str(Nt) ' terms'],'f_o(x)');
xlabel('x')
ylabel('y')

pause(0.1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End Sheet 2 Q3(a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sheet 2 Q3(b)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% Enter coefficients of cubic f(x) = Ax^3+Bx^2+Cx+D for 0 <= x <= 1
A = 1;
B = 1;
C = 1;
D = 1;

% truncate after 2*Nt+1 terms 
Nt = 5;

% plot domain
a = -3;
b = 3;

% symbolic variables
syms x n

% where n is a natural number
assume(n,'integer');
assumeAlso(n > 0);

% f_e and f_o as symbolic functions for |x| <= 1
syms fe(x) fo(x)
fe(x) = A*sign(x)*x^3 + B*x^2 + C*sign(x)*x + D;
fo(x) = A*x^3 + B*sign(x)*x^2 + C*x + D*sign(x);

% evaluate Fourier coefficients
disp('Fourier coefficients');
fe_a0 = simplify( int(fe,-1,1) )
fe_an = simplify( int(fe*cos(n*pi*x),-1,1) )
fe_bn = simplify( int(fe*sin(n*pi*x),-1,1) )
%
fo_a0 = simplify( int(fo,-1,1) )
fo_an = simplify( int(fo*cos(n*pi*x),-1,1) )
fo_bn = simplify( int(fo*sin(n*pi*x),-1,1) )

% construct fourier series as symbolic expressions

Sfe = fe_a0/2;
Sfo = fo_a0/2;

for k=1:Nt
    Sfe = Sfe + subs(fe_an,n,k)*cos(k*pi*x) + subs(fe_bn,n,k)*sin(k*pi*x);
    Sfo = Sfo + subs(fo_an,n,k)*cos(k*pi*x) + subs(fo_bn,n,k)*sin(k*pi*x);
end

% turn it into a matlab function for plotting
mSfe = matlabFunction(Sfe);
mSfo = matlabFunction(Sfo);

% sample points in domain
xx = linspace(a,b,1000);

% sample values of Fourier series
valSfe = mSfe(xx);
valSfo = mSfo(xx);

% sample values of f
% mod operation to sample values outside interval -pi < x <= pi
valfe = fe( mod(xx-1,2)-1 );
valfo = fo( mod(xx-1,2)-1 );

% plot
figure(2);
set(gcf,'Position',[0 0 900 600],'Color','w');
%
subplot(2,1,1)
plot(xx,valSfe,'r-','LineWidth',1.25);
hold on;
plot(xx,valfe,'k');
grid on;
legend(['Cosine series, ' num2str(Nt+1) ' terms'],'f_e(x)');
xlabel('x')
ylabel('y')
title(['f(x) = Ax^3+Bx^2+Cx+D for 0 <= x <= 1 with '...
        ' A = ' num2str(A)...
        ', B = ' num2str(B)...
        ', C = ' num2str(C)...
        ', D = ' num2str(D)]);
%
subplot(2,1,2)
plot(xx,valSfo,'r-','LineWidth',1.25);
hold on;
plot(xx,valfo,'k');
grid on;
legend(['Sine series, ' num2str(Nt) ' terms'],'f_o(x)');
xlabel('x')
ylabel('y')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End Sheet 2 Q3(b)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%