89 lines
2 KiB
Mathematica
89 lines
2 KiB
Mathematica
|
%% Assignment 2
|
||
|
|
% Name: Claudio Maggioni
|
||
|
|
%
|
||
|
|
% Date: 19/3/2019
|
||
|
|
%
|
||
|
|
% This is a template file for the first assignment to get started with running
|
||
|
|
% and publishing code in Matlab. Each problem has its own section (delineated
|
||
|
|
% by |%%|) and can be run in isolation by clicking into the particular section
|
||
|
|
% and pressing |Ctrl| + |Enter| (evaluate current section).
|
||
|
|
%
|
||
|
|
% To generate a pdf for submission in your current directory, use the following
|
||
|
|
% three lines of code at the command window:
|
||
|
|
%
|
||
|
|
% >> options.format = 'pdf'; options.outputDir = pwd; publish('assignment3.m', options)
|
||
|
|
%
|
||
|
|
|
||
|
|
%% 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);
|
||
|
|
|
||
|
|
%% Problem 6
|
||
|
|
[X,Y] = meshgrid((0:2000)/2000);
|
||
|
|
A = exp(-sqrt((X-Y).^2));
|
||
|
|
L = outerProductCholesky(A);
|
||
|
|
disp(norm(L*L'-A, 'fro'));
|
||
|
|
|
||
|
|
%% 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;
|
||
|
|
end
|
||
|
|
|
||
|
|
%% Problem 6 (continued)
|
||
|
|
|
||
|
|
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));
|
||
|
|
end
|
||
|
|
end
|