hw4: ex1 and ex3 done

hw4/assignment4.m

@@ -1,4 +1,4 @@
-%% Assignment 2
+%% Assignment 4
 % Name: Claudio Maggioni
 %
 % Date: 19/3/2019
@@ -15,74 +15,76 @@
 %
 
 %% Problem 3
-format long
-A = [4 3 2 1; 8 8 5 2; 16 12 10 5; 32 24 20 11];
-[L,U,P] = pivotedOuterProductLU(A);
-display(L);
-display(U);
-display(P);
+clear;
 
-%% Problem 6
-[X,Y] = meshgrid((0:2000)/2000); 
-A = exp(-sqrt((X-Y).^2));
-L = outerProductCholesky(A);
-disp(norm(L*L'-A, 'fro'));
+n=10;
+f = @(x) (exp(-(x^2)/2)/sqrt(2*pi));
+x_5 = computeEquidistantXs(5);
+x_10 = computeEquidistantXs(10);
+x_5c = computeChebyshevXs(5);
+x_10c = computeChebyshevXs(10);
 
-%% Problem 3 (continued)
-
-function [L,U,P] = pivotedOuterProductLU(A)
-    dimensions = size(A);
-    n = dimensions(1);
-    
-    p = 1:n;
-    L = zeros(n);
-    U = zeros(n);
-
-    for i = 1:n
-        values = A(:,i);
-        values(values == 0) = -Inf;
-        [~, p_k] = max(values);
-        k = find(p == p_k);
-        
-        if k == -Inf
-            disp("Matrix is singular");
-            L = [];
-            U = [];
-            P = [];
-            return
-        end
-        
-        p(k) = p(i);
-        p(i) = p_k;
-        
-        L(:,i) = A(:,i) / A(p(i),i);
-        U(i,:) = A(p(i),:);
-        A = A - L(:,i) * U(i,:);
-    end
-    
-    I = eye(n);
-    P = zeros(n);
-    for i = 1:n
-        P(:,i) = I(:,p(i));
-    end
-    P = transpose(P);
-    L = P * L;   
+xe = (-1:0.01:1)';
+y = zeros(size(xe));
+for i = 1:size(xe, 1)
+    y(i) = f(xe(i));
 end
 
-%% Problem 6 (continued)
+p_5 = computePolypoints(f, xe, x_5, 5);
+p_10 = computePolypoints(f, xe, x_10, 10);
 
-function [L] = outerProductCholesky(A)
-    if ~all(eig(A) >= 0)
-        disp("matrix is not positive definite");
-        L = [];
-        return;
-    end
-    dimensions = size(A);
-    n = dimensions(1);
-    L = zeros(n);
-    
-    for i = 1:n
-        L(:, i) = A(:, i) / sqrt(A(i,i));
-        A = A - L(:, i) * transpose(L(:, i));
+p_5c = computePolypoints(f, xe, x_5c, 5);
+p_10c = computePolypoints(f, xe, x_10c, 10);
+
+figure;
+subplot(1,2,1);
+plot(xe, p_5, xe, p_10, xe, y);
+subplot(1,2,2);
+plot(xe, p_5c, xe, p_10c, xe, y);
+
+function x = computeEquidistantXs(n)
+    x = zeros(2*n+1,1);
+    for i = 1:2*n+1
+        x(i) = (i-n-1)/n;
+    end
+end
+
+
+function x = computeChebyshevXs(n)
+    x = zeros(2*n+1,1);
+    for i = 1:2*n+1
+        x(i) = cos((2*i - 1) *pi / (4*n + 2));
+    end
+end
+
+
+function p = computePolypoints(f, xe, x, n)
+    p = zeros(size(xe));
+
+    for i = 1:(2*n+1)
+        e_i = zeros(2*n+1, 1);
+        e_i(i) = 1;
+        N = NewtonInterpolation(x, e_i);
+        p = p + f(x(i)) * HornerNewton(N, x, xe);
+    end
+end
+
+% Assuming x and y are column vectors with the same length
+function N = NewtonInterpolation (x,y)
+    n = size(x, 1);
+    N = y;
+    for i = 1:n
+        N(n:-1:i+1) = (N(n:-1:i+1) - N(n-1:-1:i)) ./ (x(n:-1:i+1) - x(n-i:-1:1));
+    end
+end
+
+% N is the array of coefficients
+% 
+% xi evaluation points
+function p = HornerNewton(N, x, xe)
+    n = size(x, 1);
+    p = ones(size(xe, 1), 1) * N(n);
+    for i = n-1:-1:1 
+        p = p .* (xe - x(i)) + N(i);
     end
 end