function ch07_oscillatory
%==========================================================================
% ch07_oscillatory.m
%
% Numerical integration of a rapidly oscillating function.
%
% The program evaluates
%
%          1
%          /
%         |  exp(-x) cos(omega*x) dx,       omega = 50,
%         /
%         0
%
% using:
%   1. Composite Simpson rule
%   2. Composite quadratic Filon rule
%
% Simpson resolves the complete oscillatory integrand numerically.
% Filon approximates only the slowly varying amplitude by a quadratic
% polynomial and integrates the oscillatory factor analytically.
%
% Based on the Fortran implementation by Alexander Godunov.
% MATLAB version prepared for the companion website, 2026.
%
% A plot of the integrand (or integration domain) is shown before the
% numerical calculation to make the character of the problem visible.
%==========================================================================

a = 0.0;
b = 1.0;
omega = 50.0;

exact = ( exp(-b)*(-cos(omega*b) + omega*sin(omega*b)) ...
        - exp(-a)*(-cos(omega*a) + omega*sin(omega*a)) ) ...
        /(1.0 + omega^2);

n_simpson = [100, 200, 400, 800, 1600];
n_filon = [20, 40, 80];

% === Plot the integrand =================================================
xplot = linspace(a,b,4001);
amplitude = f(xplot);
oscillatory = amplitude.*cos(omega*xplot);

figure('Name','Rapidly oscillating integrand')
plot(xplot,oscillatory,'LineWidth',1.2)
hold on
plot(xplot, amplitude,'--','LineWidth',1.0)
plot(xplot,-amplitude,'--','LineWidth',1.0)
hold off
xlabel('x')
ylabel('f(x) cos(\omega x)')
title('Rapidly oscillating integrand and amplitude envelope')
legend('exp(-x) cos(50x)','+ exp(-x)','- exp(-x)', ...
       'Location','best')
grid on
% === End plot ===========================================================

fprintf('Rapidly oscillating integral\n');
fprintf('Integral of exp(-x) cos(omega*x), omega = %6.1f\n', omega);
fprintf('Exact value          = %15.7e\n\n', exact);

fprintf('Composite Simpson rule\n');
fprintf('        n            Integral       Absolute error\n');
fprintf('----------------------------------------------------\n');

for n = n_simpson
    integral = composite_simpson_oscillatory(@f, a, b, omega, n);
    error_value = abs(integral-exact);
    fprintf('%9d%20.7e%20.7e\n', n, integral, error_value);
end

fprintf('\nComposite quadratic Filon rule\n');
fprintf('        n            Integral       Absolute error\n');
fprintf('----------------------------------------------------\n');

for n = n_filon
    integral = filon_quadratic(@f, a, b, omega, n);
    error_value = abs(integral-exact);
    fprintf('%9d%20.7e%20.7e\n', n, integral, error_value);
end

end


function y = f(x)
%--------------------------------------------------------------------------
% Slowly varying amplitude in integral f(x) cos(omega*x) dx.
%--------------------------------------------------------------------------
y = exp(-x);
end


function integral = composite_simpson_oscillatory(fun, a, b, omega, n)
%--------------------------------------------------------------------------
% Composite Simpson rule for f(x) cos(omega*x).
%--------------------------------------------------------------------------

if n < 2 || mod(n,2) ~= 0
    error('composite_simpson_oscillatory: n must be even.');
end

h = (b-a)/n;

sum_odd = 0.0;
for i = 1:2:n-1
    x = a+i*h;
    sum_odd = sum_odd + fun(x)*cos(omega*x);
end

sum_even = 0.0;
for i = 2:2:n-2
    x = a+i*h;
    sum_even = sum_even + fun(x)*cos(omega*x);
end

integral = h*(fun(a)*cos(omega*a) + fun(b)*cos(omega*b) ...
           + 4.0*sum_odd + 2.0*sum_even)/3.0;

end


function integral = filon_quadratic(fun, a, b, omega, n)
%--------------------------------------------------------------------------
% Composite quadratic Filon rule.
%
% On each pair of subintervals, the slowly varying amplitude is
% approximated by a quadratic polynomial. The polynomial is integrated
% analytically against cos(omega*x).
%--------------------------------------------------------------------------

if n < 2 || mod(n,2) ~= 0
    error('filon_quadratic: n must be a positive even integer.');
end

if abs(omega) < 1.0e-12
    error('filon_quadratic: omega is too close to zero.');
end

h = (b-a)/n;
theta = omega*h;

c0 = 2.0*sin(theta)/omega;

c2 = 2.0*h^2*sin(theta)/omega ...
   + 4.0*h*cos(theta)/(omega^2) ...
   - 4.0*sin(theta)/(omega^3);

s1 = -2.0*h*cos(theta)/omega ...
   +  2.0*sin(theta)/(omega^2);

integral = 0.0;

for k = 0:2:n-2

    x0 = a+k*h;
    x1 = x0+h;
    x2 = x0+2.0*h;

    f0 = fun(x0);
    f1 = fun(x1);
    f2 = fun(x2);

    a0 = f1;
    a1 = (f2-f0)/(2.0*h);
    a2 = (f0-2.0*f1+f2)/(2.0*h^2);

    pair_integral = cos(omega*x1)*(a0*c0 + a2*c2) ...
                  - sin(omega*x1)*a1*s1;

    integral = integral + pair_integral;
end

end