%=======================================================================
% Chapter 10, code 2: a linear BVP with different boundary conditions.
% Solve y'' + p(x)*y' + q(x)*y = r(x) on a uniform grid.
% Method: centered finite differences and the Thomas algorithm.
% Derivative boundaries use ghost-point elimination. All endpoints are
% included as unknowns in the tridiagonal system.
%
% Example: dimensionless steady heat conduction in a rod,
%   T'' = -sin(pi*x), 0 <= x <= 1.
% Default: T(0) = 0, T'(1) = 0 (an insulated right end).
% y denotes T. Exact solutions are supplied for all four boundary keys.
%
% MATLAB alternative: assemble a sparse matrix A and solve with y = A\rhs.
% The boundary rows must implement the selected conditions below.
%
% To solve another problem, edit the inputs and coefficients;
% replace exact_solution or set check_exact = false.
%=======================================================================
clear;

%-------------------- Problem and grid: edit here -----------------------
a = 0.0;
b = 1.0;
n = 20;                      % Number of grid intervals
key = 2;
left_value = 0.0;
right_value = 0.0;
robin_alpha = 1.0;
robin_beta = 1.0;
check_exact = true;
make_plot = true;
%
% Boundary key and meaning of the two boundary data:
%   1: y(a)  = left_value; y(b)  = right_value
%   2: y(a)  = left_value; y'(b) = right_value
%   3: y'(a) = left_value; y(b)  = right_value
%   4: y(a)  = left_value;
%      robin_alpha*y(b) + robin_beta*y'(b) = right_value
%
% Derivatives mean dy/dx, not the outward normal derivative.
% For a rod's convecting right end, key 4 expresses
%   k*T'(b) + H*T(b) = H*T_environment.
% Use dimensionless conductivity k and heat-transfer coefficient H here.
%-----------------------------------------------------------------------
if n < 2 || n ~= floor(n) || b <= a
    error('Use an integer n >= 2 and b > a.');
end
if ~ismember(key,1:4)
    error('Choose a boundary key from 1 to 4.');
end
if key == 4
    if robin_beta == 0
        error('Use key 1 for a fixed right-end value.');
    end
    if check_exact && abs(robin_alpha*(b-a)+robin_beta) <= ...
            eps*max(abs(robin_alpha*(b-a)),abs(robin_beta))
        error('The supplied example has no unique solution for these Robin data.');
    end
end

h = (b-a)/n;
x = linspace(a,b,n+1)';
y = finite_difference(@coefficients,n,x,key,left_value,right_value, ...
                      robin_alpha,robin_beta);

filename = 'ch10_boundary_conditions.dat';
fid = fopen(filename,'w');
if fid == -1
    error('Cannot open the data file.');
end

fprintf('Linear boundary value problem: boundary-condition choices\n');
fprintf('Boundary key:   %8d\n',key);
fprintf('Grid intervals: %8d\n',n);
fprintf('Grid spacing:   %12.6f\n',h);
switch key
    case 1
        fprintf('Fixed values at both ends.\n');
    case 2
        fprintf('Fixed left value; prescribed right derivative.\n');
    case 3
        fprintf('Prescribed left derivative; fixed right value.\n');
    case 4
        fprintf('Fixed left value; right Robin condition.\n');
end
fprintf('\n');

if check_exact
    exact = exact_solution(x,a,b,key,left_value,right_value,robin_alpha,robin_beta);
    abs_error = abs(y-exact);
    fprintf('%10s %16s %16s %12s\n', ...
            'x','numerical y','exact y','abs error');
    fprintf(fid,'# x numerical_y exact_y abs_error\n');
    for i = 1:n+1
        fprintf('%10.6f %16.8e %16.8e %12.4e\n', ...
                x(i),y(i),exact(i),abs_error(i));
        fprintf(fid,' %.16e %.16e %.16e %.16e\n', ...
                x(i),y(i),exact(i),abs_error(i));
    end
else
    fprintf('%10s %16s\n','x','numerical y');
    fprintf(fid,'# x numerical_y\n');
    for i = 1:n+1
        fprintf('%10.6f %16.8e\n',x(i),y(i));
        fprintf(fid,' %.16e %.16e\n',x(i),y(i));
    end
end
fclose(fid);

fprintf('\n');
if check_exact
    fprintf('Largest grid error: %12.4e\n',max(abs_error));
end
fprintf('Data file: %s\n',filename);

if make_plot
    figure;
    if check_exact
        subplot(2,1,1);
    end
    plot(x,y,'ko-','LineWidth',1.0,'MarkerSize',4);
    hold on;
    if check_exact
        x_plot = linspace(a,b,401)';
        plot(x_plot,exact_solution(x_plot,a,b,key,left_value,right_value,robin_alpha,robin_beta), ...
             'k--','LineWidth',1.2);
        legend('Finite differences','Exact solution','Location','best');
    end
    hold off;
    xlabel('x');
    ylabel('y(x)');
    title('Linear boundary value problem');
    grid on;

    if check_exact
        subplot(2,1,2);
        plot(x,abs_error,'ko-','LineWidth',1.0,'MarkerSize',4);
        xlabel('x');
        ylabel('Absolute error');
        grid on;
    end
end


function [p,q,r] = coefficients(x)
% Supply y'' + p(x)*y' + q(x)*y = r(x).
% Here y is dimensionless temperature in a uniformly conducting rod.
p = 0.0;
q = 0.0;
r = -sin(pi*x);
end


function y = exact_solution(x,a,b,key,left_value,right_value,alpha,beta)
% For y'' = -sin(pi*x), y = sin(pi*x)/pi^2 + c*(x-a) + d.
% Determine c and d from the selected boundary conditions.
% Default key 2 on [0,1]: y = sin(pi*x)/pi^2 + x/pi.
pa = sin(pi*a)/pi^2;
pb = sin(pi*b)/pi^2;
switch key
    case 1
        d = left_value-pa;
        c = (right_value-pb-d)/(b-a);
    case 2
        d = left_value-pa;
        c = right_value-cos(pi*b)/pi;
    case 3
        c = left_value-cos(pi*a)/pi;
        d = right_value-pb-c*(b-a);
    case 4
        d = left_value-pa;
        c = (right_value-alpha*(pb+d)-beta*cos(pi*b)/pi) ...
            /(alpha*(b-a)+beta);
end
y = sin(pi*x)/pi^2+c*(x-a)+d;
end


function y = finite_difference(fcn,n,x,key,left_value,right_value,alpha,beta)
%=======================================================================
% Interior equations are identical to Code 1. All n+1 values, including
% endpoints, are unknowns here. MATLAB index 1 is the left endpoint;
% index n+1 is the right endpoint.
%
% A fixed endpoint gets a row containing y = prescribed value.
% For a derivative boundary, introduce a ghost point outside the grid:
%   y_ghost_left  = y_next - 2*h*y'(a)
%   y_ghost_right = y_previous + 2*h*y'(b).
% Substitute into the centered ODE at the endpoint to remove the ghost.
% For a right Robin condition, use y'(b) = (right_value-alpha*y(b))/beta.
%=======================================================================
h = (x(end)-x(1))/n;
m = n+1;
lower = zeros(m,1);
diagonal = zeros(m,1);
upper = zeros(m,1);
rhs = zeros(m,1);

for i = 2:n
    [p,q,r] = fcn(x(i));
    lower(i) = 1-h*p/2;
    diagonal(i) = -2+h*h*q;
    upper(i) = 1+h*p/2;
    rhs(i) = h*h*r;
end

% Left boundary: derivative for key 3; fixed value otherwise.
if key == 3
    [p,q,r] = fcn(x(1));
    diagonal(1) = -2+h*h*q;
    upper(1) = 2;
    rhs(1) = h*h*r+(2*h-h*h*p)*left_value;
else
    diagonal(1) = 1;
    rhs(1) = left_value;
end

% Right boundary: fixed value, derivative, or Robin condition.
switch key
    case {1,3}
        diagonal(m) = 1;
        rhs(m) = right_value;
    case 2
        [p,q,r] = fcn(x(m));
        lower(m) = 2;
        diagonal(m) = -2+h*h*q;
        rhs(m) = h*h*r-(2*h+h*h*p)*right_value;
    case 4
        [p,q,r] = fcn(x(m));
        factor = (2*h+h*h*p)/beta;
        lower(m) = 2;
        diagonal(m) = -2+h*h*q-factor*alpha;
        rhs(m) = h*h*r-factor*right_value;
end

y = thomas(lower,diagonal,upper,rhs);
end


function y = thomas(lower,diagonal,upper,rhs)
%=======================================================================
% Solve a tridiagonal system by forward elimination and backward
% substitution. The local copies of diagonal and rhs are overwritten.
%
% No pivoting is used. A near-zero pivot stops the calculation.
% Such a failure can require a solver with pivoting or can indicate
% a singular boundary value problem.
%=======================================================================
m = length(rhs);
y = zeros(m,1);

% Check each pivot before division.
for i = 1:m
    scale = max([abs(lower(i)),abs(diagonal(i)),abs(upper(i))]);
    if abs(diagonal(i)) <= eps*scale
        error('Thomas algorithm: zero or near-zero pivot.');
    end
    if i == m
        break;
    end
    multiplier = lower(i+1)/diagonal(i);
    diagonal(i+1) = diagonal(i+1)-multiplier*upper(i);
    rhs(i+1) = rhs(i+1)-multiplier*rhs(i);
end

y(m) = rhs(m)/diagonal(m);
for i = m-1:-1:1
    y(i) = (rhs(i)-upper(i)*y(i+1))/diagonal(i);
end
end