.

Sheet6/Ex6A.m

@@ -0,0 +1,48 @@
+clear; clc; close all;
+pList = [5, 10, 20, 100];
+nIter = 9;
+
+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);
+
+        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')