Sheet 4 - solutions

Sheet_4/Ex_4_B.m

@@ -0,0 +1,32 @@
+% 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');

Sheet_4/Ex_4_C.m

@@ -0,0 +1,51 @@
+% 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');