%% Practical 1 -- Model solutions
% Asgeir Birkisson and Mohsin Javed

clear, clc, close all

%% Problem 2
% Create a vector to be put on the sub- and super diagonals
e = ones(15, 1);

%%
% Call the diag method with the vector e, and the location we want it on (-1 for
% subdiagonal, +1 for superdiagonal, no extra argument for diagonal). Notice for
% the diagonal, we need to add an entry to e to get the correct dimensions.
D2 = diag(e, -1) + diag(-2*[e; 1]) + diag(e, 1);

%%
% Manually modify the top-right and bottom-left entries.
D2(1,16) = 1; 
D2(16,1) = 1;
% Scale the matrix
h = 2/16;
D2 = (1/h)^2 * D2;

%% Problem 3
% Do the same as above, using the >>toeplitz<< method.
D2 = (1/h)^2 * toeplitz([-2 1 zeros(1, 13) 1]);

%% Problem 4a
% Create variables for N and h that we can easily adjust
N = 16;
h = 2/N;

%%
% Do same as in Problem 2, but now using N to control the size of the vector e
% and the matrices.
e = ones(N-1, 1);
D2 = diag(e, -1) + diag(-2*[e; 1]) + diag(e, 1);
D2(1, N) = 1; 
D2(N, 1) = 1;
D2 = (1/h)^2 * D2;

%% Problem 4b
% Same as above, but now using the toeplitz method
D2 = (1/h)^2 * toeplitz([-2 1 zeros(1, N-3) 1]);

%% Problem 5
% Define a vector x, and use the sin command on x to define f
x = -1:2/N:1-2/N;

%%
% Transpose x to give a column vector
x = x';
f = sin(pi*x);
D2f = D2 * f;
plot(x, D2f)

%% Problem 6
% Plot the finite difference approximation to the second derivative of the 
% function sin(pi*x).
figure(1)
subplot(2, 1, 1)
plot(x, D2f);
title('Finite difference approximation to d^2/dx^2 sin(\pi x)');
xlabel('x');
ylabel('-\pi^2 sin(\pi x)');

%%
% Add a plot showing the error
subplot(2, 1, 2)
plot(x, (-pi^2*f) - D2f);
xlabel('x')
ylabel('pointwise error');

%% Problem 7
clf;
N = 256;
h = 2/N;
D2 = (1/h)^2 * toeplitz([-2 1 zeros(1, N-3) 1]);
x = (-1:2/N:1-2/N).';
f = exp(sin(pi*x));
fpp = @(x) pi^2 * exp(sin(pi*x)) .* cos(pi*x).^2 ...
         - pi^2 * exp(sin(pi*x)) .* sin(pi*x)

plot(x, f)
hold on
plot(x, D2*f)
plot(x, fpp(x))
hold off
%%
figure
plot(x, D2*f - fpp(x))

%% Problem 8
% Set up a second derivative finite difference matrix
N = 256; 
h = 2/N;
D2 = (1/h)^2 * toeplitz([-2 1 zeros(1, N-3) 1]);

%%
% Subtract a finite difference approximation of the identity operator, and scale
% the matrix.
A = (D2-eye(N))/(1+pi^2);

%%
% Set up a right-hand side
x = -1 : 2/N : 1-2/N;
x = x';
f = sin(pi*x);

%%
% Solve the linear system
y = A\f;

% Plot the results
figure(3)
subplot(2, 1, 1)
plot(x, y);
title('y'''' - y = (1+\pi^2)*sin(\pi x), y(-1) = y(1), y''(-1) = y''(1) ');
xlabel('x');
ylabel('Solution: y(x) = -sin(\pi x)');

subplot(2, 1, 2)
plot(x, y - (-sin(pi*x)));
xlabel('x')
ylabel('pointwise error');

%% Problem 9
% Same as above, but using sparse matrices
N = 256;
h = 2/N;
e = ones(N, 1);
D2 = spdiags([e -2*e e], -1:1, N, N);
D2(1,N) = 1;
D2(N,1) = 1;
D2 = (1/h)^2 * D2;

x = -1:2/N:1-2/N;
x = x';
f = sin(pi*x);
D2f = D2*f;

figure(4)
subplot(2, 1, 1)
plot(x, D2f);
title('Finite difference approximation to d^2/dx^2 sin(\pi x)');
xlabel('x');
ylabel('-\pi^2 sin(\pi x)');

subplot(2, 1, 2)
plot(x, (-pi^2*f) - D2f);
xlabel('x')
ylabel('pointwise error');

%% Problem 10
% Find eigenvalues of the matrix D2
N = 2^6;
h = 2/N;
e = ones(N, 1);
D2 = spdiags([e -2*e e], -1:1, N, N);
D2(1,N) = 1;
D2(N,1) = 1;
D2 = (1/h)^2 * D2;
eigs(D2)
plot(eig(D2),'.','MarkerSize',10)
n=N/2;
ev=2*(cos((0:n)*pi*h)-1)/h^2;
ev=[ev, ev(2:end-1)];
ev=sort(ev);
hold on
plot(ev,'rx','MarkerSize',10)
figure
plot(ev(:)-eig(D2))