Uploading ex6

Sheet_6/adapt_h.m

@@ -0,0 +1,18 @@
+% h-adaptivity
+% for a given mesh and solution we generate a new adapted mesh by splitting
+% the elements with large errors into two new elements
+
+function mesh = adapt_h(nodes,u,lambda)
+    n_nodes = length(nodes);
+    flux_jumps = jumps_flux(nodes,u,lambda);
+    alpha = 0.5;    % parameter for choosing elements to refine
+    crit_error = alpha*max(abs(flux_jumps));
+    mesh = nodes;   % will be the new nodes/mesh
+    % finding elements to refine
+    for i=2:n_nodes-1
+        if abs(flux_jumps(i)) > crit_error
+            mesh = [mesh, nodes(i-1)+(nodes(i)-nodes(i-1))/2, nodes(i)+(nodes(i+1)-nodes(i))/2];
+        end
+    end
+    mesh = unique(mesh);
+end

Sheet_6/adapt_r.m

@@ -0,0 +1,24 @@
+% r-adaptivity
+% for a given mesh and solution we generate a new adapted mesh by moving
+% the existing nodes within the mesh; they get moved to a position in order
+% to equally distribute the error over the intervall
+% we use the De Boor's algorithm (Huang, Russell; Adaptive Moving Mesh
+% Methods; $ 2.2.1)
+
+function mesh = adapt_r(nodes,u,lambda)
+    n_nodes = length(nodes);
+    flux_jumps = abs(jumps_flux(nodes,u,lambda));
+    p = 1/2*(flux_jumps(1:end-1) + flux_jumps(2:end));  % has values for each element
+    h_vec = nodes(2:end) - nodes(1:end-1);
+    P = cumsum(h_vec'.*p);
+    P = [0;P];
+    xi = linspace(0,1,n_nodes);
+    mesh = nodes;
+    for j=2:n_nodes-1
+        idx_k = find(xi(j)*P(end) <= P, 1, 'first');
+        if idx_k > 1
+            x_j = nodes(idx_k-1) + xi(j)*(P(end)-P(idx_k -1))/p(idx_k -1);
+            mesh(j) = x_j;
+        end
+    end
+end

Sheet_6/assembling.m

@@ -0,0 +1,20 @@
+% assembling of the stiffness matrix and the load vector for a given mesh
+% in 1D
+
+function [K,f_vec] = assembling(nodes,lambda,f)
+    n_nodes = length(nodes);
+    K = zeros(n_nodes);
+    f_vec = zeros(n_nodes,1);
+    h_diff_vec = nodes(2:end) - nodes(1:end-1); % step width
+    for k = 1:n_nodes-1
+        % stiffness matrix
+        lambda_int = integral(lambda, nodes(k), nodes(k+1));
+        K_loc = lambda_int/h_diff_vec(k)^2*[1,-1;-1,1];
+        K(k:k+1,k:k+1) = K(k:k+1,k:k+1) + K_loc;
+        % right hand side
+        phi_left = @(x) (nodes(k+1)-x)/h_diff_vec(k).*f(x);
+        phi_right = @(x) (x-nodes(k))/h_diff_vec(k).*f(x);
+        f_loc = [integral(phi_left,nodes(k),nodes(k+1)); integral(phi_right,nodes(k),nodes(k+1))];
+        f_vec(k:k+1) = f_vec(k:k+1) + f_loc;
+    end
+end

Sheet_6/jumps_flux.m

@@ -0,0 +1,16 @@
+% calculation of the flux jumps according to the chosen mesh
+% nodes - corresponding to the mesh
+% u - approximation
+% flux_jumps - jump of the flux in each node of the mesh including the
+% paramter function lambda
+
+function flux_jumps = jumps_flux(nodes, u, lambda)
+    n_nodes = length(nodes);
+    flux_jumps = zeros(n_nodes,1);
+    diff_u = u(2:end) - u(1:end-1);
+    diff_x = nodes(2:end) - nodes(1:end-1);
+    eps = 1e-8;
+    for i = 2:n_nodes-1
+        flux_jumps(i) = diff_u(i)/diff_x(i)*lambda(nodes(i)+eps) - diff_u(i-1)/diff_x(i-1)*lambda(nodes(i)-eps);
+    end
+end