50 lines
1.2 KiB
Matlab
50 lines
1.2 KiB
Matlab
%% FEM vs analytical solution for -u'' + a*u = f, mixed BC
|
|
clear; clc;
|
|
%% Parameters
|
|
n = 10; % number of elements
|
|
a = 2; % reaction term
|
|
alpha = 5; % Robin BC
|
|
gB = 1; % Robin BC value
|
|
f = 3; % constant RHS
|
|
h = 1/n;
|
|
x = linspace(0,1,n+1);
|
|
|
|
%% Assemble FEM matrices
|
|
A = zeros(n+1,n+1);
|
|
F = zeros(n+1,1);
|
|
|
|
% Linear shape functions, 1D FEM
|
|
K_diff = (1/h)*[1,-1;-1,1];
|
|
K_reac = a*(h/6)*[2,1;1,2];
|
|
|
|
for e=1:n
|
|
nodes = [e, e+1];
|
|
A(nodes,nodes) = A(nodes,nodes) + K_diff + K_reac;
|
|
F(nodes) = F(nodes) + f*(h/2)*[1;1]; % f*v integral approx
|
|
end
|
|
|
|
% Dirichlet BC u(0)=0
|
|
A(1,:) = 0; A(:,1) = 0; A(1,1) = 1;
|
|
F(1) = 0;
|
|
|
|
% Robin BC at x=1
|
|
A(end,end) = A(end,end) + alpha;
|
|
F(end) = F(end) + alpha*gB;
|
|
|
|
%% Solve FEM
|
|
u_fem = A\F;
|
|
|
|
%% Analytical solution
|
|
syms C
|
|
C1 = (alpha*gB - (f/a)*alpha*(1 - exp(-sqrt(a))) - (f/sqrt(a))*exp(-sqrt(a))) ...
|
|
/ (alpha*sinh(sqrt(a)) + sqrt(a)*cosh(sqrt(a)));
|
|
u_exact = @(x) C1*sinh(sqrt(a)*x) + (f/a)*(1 - exp(-sqrt(a)*x));
|
|
u_an = u_exact(x);
|
|
|
|
%% Plot
|
|
plot(x,u_fem,'o-','LineWidth',2); hold on;
|
|
plot(x,u_an,'r--','LineWidth',2);
|
|
legend('FEM','Analytical');
|
|
xlabel('x'); ylabel('u(x)');
|
|
title('FEM vs Analytical solution');
|
|
grid on;
|