52 lines
1.3 KiB
Matlab
52 lines
1.3 KiB
Matlab
clear; clc; close all;
|
|
pList = [5, 10, 20, 100];
|
|
nIter = 5;
|
|
|
|
figure;
|
|
tiledlayout(2,2)
|
|
for ip = 1:length(pList)
|
|
p = pList(ip);
|
|
|
|
f = @(x) 2*p^3*x./(p^2*x.^2+1).^2; %RHS
|
|
lambda = @(x) 1; %DIffusion coefficient
|
|
|
|
u_exact = @(x) atan(p*x);
|
|
|
|
x = linspace(-1,1,10); %initial uniform mesh (includes x=0)
|
|
|
|
for it = 1:nIter %adaptive FEM loop
|
|
[K,F] = assemble_1D(x,lambda,f); %assemble global stiffness matrix K and load vector F
|
|
%Neumann boundary conditions at x=1
|
|
F(end) = F(end) + p/(p^2+1);
|
|
%Dirichlet boundary conditions at x=-1
|
|
K(1,:) = 0;
|
|
K(1,1) = 1;
|
|
F(1) = -atan(p);
|
|
% GH: Dirichlet boundary conditions at x=+1
|
|
K(end,:) = 0;
|
|
K(end,end) = 1;
|
|
F(end) = atan(p);
|
|
|
|
u = K \ F;
|
|
|
|
if it < nIter % h-adaptivity loop
|
|
eta = flux_jump(x,u,lambda);
|
|
x = h_adapt(x,eta,0.3);
|
|
end
|
|
end
|
|
|
|
u_ex = u_exact(x(:)); % force column vector
|
|
err_inf = max(abs(u - u_ex));
|
|
nexttile
|
|
plot(x,u,'-o','LineWidth',1.5); hold on
|
|
plot(x,u_ex,'--','LineWidth',1.5)
|
|
title(['p = ',num2str(p)])
|
|
xlabel('x'), ylabel('u')
|
|
legend('FEM','exact','Location','best')
|
|
grid on
|
|
|
|
fprintf('p = %d\n', p);
|
|
fprintf('error = %.3e\n\n', err_inf);
|
|
end
|
|
|
|
sgtitle('Exercise 6A: h-adaptivity and exact solution comparison')
|