Numerical Analysis---Hermite Interpolation

来源：互联网发布：java 泛型 t 编辑：程序博客网时间：2024/05/16 13:43

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

昨天把数值分析的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的%%%%%%%%%%%

ticN = 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 columnstriangle = 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));endx_ = 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 intervalfor inx = 1 : 100     for k = 2 : N        x_val(inx,k) = x_val(inx,k - 1) * (x_(inx) - triangle(k-1,1));    end    endy = x_val * coeff';hold onp1 = 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);endfor i = 2 : 2 : N-2    grid(i,3) = (grid(i+1,2) - grid(i,2)) / (x(i+1) - x(i));endfor 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));    endendgridfigurehold onx_ = 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;    endendplot(x_, y, '*r');

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

开始复习。。。

0 0