%% set the parameters
a = 0;
b = 0;
N = 100; %number of interior nodes
h=1./(N+1);
x_array = 0:h:1; %define the physical domain, x = (0,1)
x_array = transpose(x_array(2:N+1 )); %remove the 0 and 1, which are the boundaries
%% construct the K matrix
K_mat = zeros(N, N);
for i=1:N
for j=1:N
if i==j
K_mat(i,j) = 2./h^2;
elseif abs(i-j)==1
K_mat(i,j) = -1/h^2;
end
end
end
%% construct the right-hand side, f
f_vec = ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY
%% solve for u
u_vec = ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY
%% Analytical solution
u_analy = ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY
figure(1)
plot(x_array,u_analy)
hold on
plot(x_array,u_vec)
xlabel("x")
ylabel("u(x)")
title("Analytical and Numerical Solution for -u''=cos(2x)")
legend("Analytical","Numerical")
%% Convergence Analysis
err_array = zeros(1,99);
for N=2:100
% N = 100; %number of interior nodes
h=1./(N+1);
x_array = 0:h:1; %define the physical domain, x = (0,1)
x_array = transpose(x_array(2:N+1 )); %remove the 0 and 1, which are the boundaries
%% construct the K matrix
K_mat = zeros(N, N);
for i=1:N
for j=1:N
if i==j
K_mat(i,j) = 2./h^2;
elseif abs(i-j)==1
K_mat(i,j) = -1/h^2;
end
end
end
%% construct the right-hand side, f
f_vec = ;% ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY
%% solve for u
u_vec = ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY
%% Analytical solution
u_analy = ;% ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY
err_array(1,N-1) = sqrt(sum((u_vec-u_analy).^2)/N); %compute the L^2 error of the solution
end
figure(2)
fit_params = polyfit(log(2:1:100),log(err_array),1); %linear fit on log-log scale
slope = fit_params(1);
y_int = fit_params(2);
loglog((2:1:100),err_array ,'o')
hold on
loglog((2:1:100), (2:1:100).^slope * exp(y_int),'r')
title("L2 Norm Decay of the Solution for -u''=cos(2x)")
xlabel("N")
ylabel("L^2 Error")
txt = ['Slope = ' num2str(slope) ];
text(10,0.0001,txt,'FontSize',14)