This repository has been archived on 2021-03-28. You can view files and clone it, but cannot push or open issues or pull requests.
ICS/hw4/assignment4.m

91 lines
2.1 KiB
Mathematica
Raw Normal View History

2020-05-14 20:23:08 +00:00
%% Assignment 4
2020-05-14 13:52:49 +00:00
% 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
2020-05-14 20:23:08 +00:00
clear;
2020-05-14 13:52:49 +00:00
2020-05-14 20:23:08 +00:00
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);
2020-05-14 13:52:49 +00:00
2020-05-14 20:23:08 +00:00
xe = (-1:0.01:1)';
y = zeros(size(xe));
for i = 1:size(xe, 1)
y(i) = f(xe(i));
end
2020-05-14 13:52:49 +00:00
2020-05-14 20:23:08 +00:00
p_5 = computePolypoints(f, xe, x_5, 5);
p_10 = computePolypoints(f, xe, x_10, 10);
2020-05-14 13:52:49 +00:00
2020-05-14 20:23:08 +00:00
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;
2020-05-14 13:52:49 +00:00
end
2020-05-14 20:23:08 +00:00
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));
2020-05-14 13:52:49 +00:00
end
end
2020-05-14 20:23:08 +00:00
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);
2020-05-14 13:52:49 +00:00
end
2020-05-14 20:23:08 +00:00
end
% Assuming x and y are column vectors with the same length
function N = NewtonInterpolation (x,y)
n = size(x, 1);
N = y;
2020-05-14 13:52:49 +00:00
for i = 1:n
2020-05-14 20:23:08 +00:00
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);
2020-05-14 13:52:49 +00:00
end
end