Numerical Analysis---Hermite Interpolation

这两天都在准备考试，没怎么动弹。

昨天把数值分析的take-home做了，第一个编程题就是Hermite Interpolation.

Hermite Interpolation需要几组点以及他们的一阶导数值，因此hermite interpolation在区间内都是连续的，这个性质很好~。例题中其结果与原函数的plot基本重合。

（x = [-1 -1 -1/2 -1/2 0 0 1/2 1/2 1 1]';

    f = @(x) 1 ./ (1 + x.^2);）

先上图感受一下~



码渣分别用loop和vector来写，原本以为vectorize之后能更快，没想到却慢了一倍的时间。。。希望大神赐教。

%%%%%%%%%用vector的%%%%%%%%%%%

tic
N = 10;
x = [-1 -1 -1/2 -1/2 0 0 1/2 1/2 1 1]';

f = @(x) 1 ./ (1 + x.^2);
df = @(x) -(2*x)./(x.^2 + 1).^2;

%initialize first three columns
triangle = zeros(N, N + 1);
triangle(:,1) = x;
triangle(:,2) = f(x);
triangle(1:2:N-1, 3) = df(x(1:2:N-1));
triangle(2:2:N-2, 3) = (triangle(3:2:N-1,2) - triangle(2:2:N-2,2)) ./ ...
(x(4:2:N) - x(2:2:N-2));
%calculate the rest of the triangle
for i = 1 : N-2
triangle(1:N-1-i,i+3) = (triangle(2:N-i,i+2) - triangle(1:N-1-i,i+2))./...
(x(i+2:N) - x(1:N-i-1));
end

x_ = linspace(-1, 1, 100);

coeff = triangle(1, 2:N+1);
x_val = ones(100, N);
y_val = zeros(1,100);

%calculat 100 more points at the interval
for inx = 1 : 100
for k = 2 : N
x_val(inx,k) = x_val(inx,k - 1) * (x_(inx) - triangle(k-1,1));
end
end
y = x_val * coeff';

hold on
p1 = plot(x_, f(x_),'-b','LineWidth',3);
p2 = plot(x_, y, '.y');
legend('f(x)','H(x)');
toc

%%%%%%%%%%%用loop的%%%%%%%%%%%%%%
N = 10;
x = [-1 -1 -1/2 -1/2 0 0 1/2 1/2 1 1]';

f = 1 ./ (1 + x.^2);
df = -(2*x)./(x.^2 + 1).^2;

grid = zeros(N, N+1);

grid(:,1) = x;
grid(:,2) = f;

for i = 1 : 2 : N-1
grid(i,3) = df(i);
end

for i = 2 : 2 : N-2
grid(i,3) = (grid(i+1,2) - grid(i,2)) / (x(i+1) - x(i));
end

for i = 1 : N-2
for j = 1 : N-1-i
grid(j,i+3) = (grid(j+1,i+2) - grid(j,i+2)) / (x(j+i+1) - x(j));
end
end

grid
figure
hold on

x_ = linspace(-1, 1, 100);
df = -(2*x_)./(x_.^2 + 1).^2;
f = 1 ./ (1 + x_.^2);
plot(x_, f)

y = grid(1,2) * ones(1,100);
temp = 1;
for inx = 1 : 100
for i = 3 : N+1
for j = 1 : i - 2
temp = temp * (x_(inx) - grid(j,1));
end
temp = temp * grid(1,i);
y(inx) = y(inx) + temp;
temp = 1;
end
end

plot(x_, y, '*r');


这学期Numerical Analysis算是最challenging的课了，准备暑假再好好总结一下，再用html5把这些算法图像化玩玩，希望自己不会食言~
开始复习。。。