插值多项式

三种等价的插值方法:多项式插值、拉格朗日插值、牛顿插值。
还有一个三次样条插值,有三种边界条件。
下面是写的matlab命令,就暂且当命令吧(一直觉得不能称之为代码)


function sth=Q2()

%{
    Rx is the Runge function
    function spline3 is the Spline interpolation function in second-order boundary condition
    function lpline is the Lagrange polynomial interpolation
    function pline is the Polynomial interpolation
    function npline is the Newton polynomial interpolation
%}


    format long;

    n = 10;

    Rx = @(x) 1./(1+25*x.^2);
    x = linspace(-1, 1, n);
    y = Rx(x);
   
    x_val = -1.0:0.001:1.0;
    S_val = spline3(x, y, 0, 0, x_val);
    L_val = lpline(x, y, x_val);
    P_val = pline(x, y, x_val);
    N_val = npline(x, y, x_val);

    fg = figure(1);

    plot(x_val, P_val, 'g', x_val, S_val, 'r', ...
        x_val, L_val, 'c--', x_val, N_val, 'b.', ...
        x, y, 'ro', 'MarkerSize', 10, 'LineWidth', 2);
    % axis([0, 1, 0, 1]);
    grid on;
    set(fg, 'Position', [0 0 1000 600]);
    title(['n = ' int2str(n)]);
    % xlabel('X');ylabel('Y');

    
    
    function S_val = spline3(x, y, f0, fn, x_val)
        format long;
        
        k = length(x);

        h = zeros(1, k-1);
        % tmp: primary difference quotient
        tmp = zeros(1, k-1);
        for i=1:(k-1)
            h(i) = (x(i+1)-x(i));
            tmp(i) = (y(i+1)-y(i))/h(i);
        end

        D = zeros(k, 1);
        for i=2:(k-1)
            D(i) = 6*(tmp(i)-tmp(i-1))/(h(i-1)+h(i));
        end
        D(1) = f0;
        D(k) = fn;

        miu = zeros(1, k-1);
        lmd = zeros(1, k-1);
        for i=1:(k-2)
            miu(i) = h(i)/(h(i)+h(i+1));
            lmd(i+1) = 1-miu(i);
        end

        T = diag(lmd, 1)+diag(miu, -1)+eye(k)*2;
        M = T\D;
        S_val = zeros(size(x_val));
        for i=1:(k-1)
            S_tmp = M(i)*(x(i+1)-x_val).^3/6/h(i)+M(i+1)*(x_val-x(i)).^3/6/h(i)...
                +(y(i)-M(i)*h(i)^2/6)*(x(i+1)-x_val)/h(i)+(y(i+1)-M(i+1)*h(i)^2/6)*(x_val-x(i))/h(i);
            if i==k-1
                S_val = S_val+S_tmp.*(x(i)<=x_val & x_val<=x(i+1));
            else
                S_val = S_val+S_tmp.*(x(i)<=x_val & x_val<x(i+1));
            end
        end
    end

    function L_val = lpline(x, y, x_val)
        k = length(x);
        d = length(x_val);
        
        x = x(:);
        y = y(:);
        x_val = x_val(:);
         
        a = ones(1, k-1);
        b = ones(d, 1);
        
        L_val = 0;
        for i=1:k
            x_tmp = x([1:i-1 i+1:k]);
            L_val = L_val+y(i)*prod((x_val*a-b*x_tmp')./(b*(x(i)*a-x_tmp')), 2);
        end
    end

    function P_val = pline(x, y, x_val)
        v = vander(x);
        a = v\y(:);
        P_val = polyval(a, x_val);
    end

    function N_val = npline(x, y, x_val)
        k = length(x);
        nc = zeros(1, k);
        nc(1) = y(1); % initial value

        for i=1:(k-1)
            tmp = zeros(1, k);
            for j=(i+1):k
                tmp(j) = (y(j)-y(j-1))/(x(j)-x(j-i));
            end
            y = tmp;
            nc(i+1) = y(i+1);
        end

        N_val = nc(1);
        for i=1:(k-1)
            N_tmp = 1;
            for j=1:i
                N_tmp = N_tmp.*(x_val-x(j));
            end
            N_val = N_val+nc(i+1)*N_tmp;
        end
    end

end

APPEND 15.11.11: 下面的程序新增了返回函数,比上面的只是返回对应于x_val的函数值要稍微显得高级一些。


function sth=Q22()
    format long;

    x = [0.2, 0.4, 0.6, 0.8, 1.0];
    y = [0.98, 0.92, 0.81, 0.64, 0.38];
    x_val = 0.2:0.001:1.0;
    
    Sx = spline3(x, y, 0, 0);
    Sx = eval(['@(xx)', vectorize(Sx)]);
    % disp(Sx);
    Lx = lpline(x, y);
    Lx = eval(['@(xx)', vectorize(Lx)]);
    %disp(Lx);
    Px = lpline(x, y);
    Px = eval(['@(xx)', vectorize(Px)]);
    %disp(Px);
    Nx = lpline(x, y);
    Nx = eval(['@(xx)', vectorize(Nx)]);
    %disp(Px);
    
    %Sx = char(Sx);
    %assignin('base','xx',1);
    %a = evalin('base',Sx);
    %disp(a);

    fg = figure(1);
    plot(x_val, Sx(x_val), 'b', x_val, Lx(x_val), 'c', ...
        x_val, Px(x_val), 'r', x_val, Nx(x_val), 'r', ...
        x, y, 'ro', 'MarkerSize', 10, 'LineWidth', 2);
    %axis([0, 1, 0, 1]);
    grid on;
    set(fg, 'Position', [0 0 1000 600]);
  
    
    function Sx = spline3(x, y, f0, fn)
        k = length(x);
        h = zeros(1, k-1);
        % tmp: primary difference quotient
        tmp = zeros(1, k-1);
        for i=1:(k-1)
            h(i) = (x(i+1)-x(i));
            tmp(i) = (y(i+1)-y(i))/h(i);
        end

        D = zeros(k, 1);
        for i=2:(k-1)
            D(i) = 6*(tmp(i)-tmp(i-1))/(h(i-1)+h(i));
        end
        D(1) = f0;
        D(k) = fn;
        % disp(D);
        
        miu = zeros(1, k-1);
        lmd = zeros(1, k-1);
        for i=1:(k-2)
            miu(i) = h(i)/(h(i)+h(i+1));
            lmd(i+1) = 1-miu(i);
        end

        T = diag(lmd, 1)+diag(miu, -1)+eye(k)*2;
        % disp(T);
        
        M = T\D;
        % disp(M);
        
        xx = 0; % inevitable step
        syms xx;
        Sx = 0;
        for i=1:(k-1)
            S_tmp = M(i)*(x(i+1)-xx).^3/6/h(i)+M(i+1)*(xx-x(i)).^3/6/h(i)...
                +(y(i)-M(i)*h(i)^2/6)*(x(i+1)-xx)/h(i)+(y(i+1)-M(i+1)*h(i)^2/6)*(xx-x(i))/h(i);
            if i==(k-1)
                Sx = Sx+S_tmp.*(x(i)<=xx & xx<=x(i+1));
            else
                Sx = Sx+S_tmp.*(x(i)<=xx & xx<x(i+1));
            end
        end
    end

    function Lx = lpline(x, y)
        k = length(x);
        x = x(:);
        y = y(:);
        a = ones(1, k-1);
        
        xx = 0; % inevitable step
        syms xx;
        Lx = 0;
        for i=1:k
            x_tmp = x([1:i-1 i+1:k]);
            Lx = Lx+y(i)*prod((xx(:)*a-ones(size(xx(:)))*x_tmp') ...
                ./(ones(size(xx(:)))*(x(i)*a-x_tmp')), 2);
        end
    end

    function Px = pline(x, y)
        v = vander(x);
        a = v\y(:);
        xx = 0; % inevitable step
        syms xx;
        Px = polyval(a, xx);
    end

    function Nx = npline(x, y)
        k = length(x);
        nc = zeros(1, k);
        nc(1) = y(1); % initial value

        for i=1:(k-1)
            tmp = zeros(1, k);
            for j=(i+1):k
                tmp(j) = (y(j)-y(j-1))/(x(j)-x(j-i));
            end
            y = tmp;
            nc(i+1) = y(i+1);
        end

        Nx = nc(1);
        xx = 0; % inevitable step
        syms xx;
        for i=1:(k-1)
            N_tmp = 1;
            for j=1:i
                N_tmp = N_tmp.*(xx-x(j));
            end
            Nx = Nx+nc(i+1)*N_tmp;
        end
    end

end