51 lines
No EOL
1.4 KiB
Matlab
51 lines
No EOL
1.4 KiB
Matlab
% Sheet 4 / Task C
|
|
clf
|
|
|
|
nvec = [10, 20, 30, 40, 70]; % vector with number of elements
|
|
p = 70; % parameter of ode
|
|
|
|
for k=1:length(nvec)
|
|
ne = nvec(k);
|
|
nodes = linspace(0,1,ne+1);
|
|
|
|
K = zeros(ne+1,ne+1); % global stiffness matrix
|
|
f = zeros(ne+1,1); % load vector
|
|
|
|
for e = 1:ne
|
|
xleft = nodes(e);
|
|
xright = nodes(e+1);
|
|
|
|
% first part of the bilinearform phi'*phi'
|
|
funcone = @(x) 1 +0*x;
|
|
Ke = ne^2*[integral(funcone,xleft,xright), -integral(funcone,xleft,xright);
|
|
-integral(funcone,xleft,xright), integral(funcone,xleft,xright)];
|
|
|
|
% second part: phi'*phi
|
|
func1 = @(x) xright - x;
|
|
func2 = @(x) x - xleft;
|
|
Ke = Ke + p*ne^2*[-integral(func1,xleft,xright), integral(func1,xleft,xright);
|
|
-integral(func2,xleft,xright), integral(func2,xleft,xright)];
|
|
|
|
K(e:e+1, e:e+1) = K(e:e+1, e:e+1) + Ke;
|
|
end
|
|
|
|
% disp(K)
|
|
|
|
% adaption for dirichlet boundary
|
|
K(1,1) = K(1,1)*(1 + 1e6);
|
|
f(end) = K(end,end)*1e6;
|
|
K(end,end) = K(end,end)*(1 + 1e6);
|
|
|
|
u = K\f;
|
|
|
|
hold on
|
|
plot(nodes,u)
|
|
title('approximation of solution u for p=70');
|
|
xlabel('x values');
|
|
ylabel('u_h');
|
|
grid on;
|
|
leg = legend(num2str(nvec'),'Location', 'Northwest');
|
|
title(leg, 'n =');
|
|
|
|
end
|
|
saveas(gcf, 'ex_4_C_plot.jpg'); |