32 lines
747 B
Matlab
32 lines
747 B
Matlab
% Sheet 4 / Task B
|
|
|
|
ne = 50; % number of elements
|
|
nodes = linspace(0,1,ne+1);
|
|
|
|
K = zeros(ne+1,ne+1); % global stiffness matrix
|
|
f = zeros(ne+1,1); % load vector
|
|
|
|
lambda = @(x) (x<1/sqrt(2))*1 + (x>=1/sqrt(2))*10;
|
|
|
|
for e = 1:ne
|
|
xleft = nodes(e);
|
|
xright = nodes(e+1);
|
|
Ke = ne^2*[integral(lambda,xleft,xright), -integral(lambda,xleft,xright);
|
|
-integral(lambda,xleft,xright), integral(lambda,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;
|
|
|
|
plot(nodes,u)
|
|
title('approximation of solution u');
|
|
xlabel('x values');
|
|
ylabel('u_h');
|