Matlab中牛顿插值看代码矩阵化的必要性

2014年的国庆假日,来的并不以往晚一些(废话)。小伙伴们都愉快的回家与家人团聚或是进城感受社会主义现代化成果。懒得出奇的我自然没有回家,也没有在这个“良辰吉日”出去玩。在宿舍宅了些许天,片子看了不少,作业倒是没动。这不假期没几天了,赶紧把作业补了。

数值分析作业里有唯一的一道代码题,很是稀疏平常,是求一个四次牛顿插值多项式。想来自己好久没写matlab了,没事动动手,要不都锈死了。

代码也挺稀疏平常,就是把书上的东西搬成了代码而已,不过原来写代码的时候没怎么注意过代码的效率问题。这次突发奇想想看看矩阵运算和循环能差多少。

先把代码贴这把。

首先是函数的m文件,NewtonPoly.m:

function [ out ] = NewtonPoly( pts, fvalue, in )
if length(pts)~=length(fvalue)
    return;
end
ptsLen = length(pts);
inLen = length(in);
bigMtx = zeros(ptsLen,ptsLen+1);
bigMtx(:,1)=pts';
bigMtx(:,2)=fvalue';
tic;
for i=1:(ptsLen-1)
    bigMtx(i+1:ptsLen,i+2)=(bigMtx(i+1:ptsLen,i+1)-bigMtx(i:ptsLen-1,i+1))./(bigMtx(i+1:ptsLen,1)-bigMtx(1:ptsLen-i,1));
end
toc;
out = zeros(1,inLen);
subX = ones(1,ptsLen);
out(1,1:inLen)=bigMtx(1,2);
tic;
for i=1:inLen
    subX(1,1:ptsLen) = in(i);
    subX = subX - pts;
    for j=2:ptsLen
        out(i)=out(i)+bigMtx(j,j+1)*prod(subX(1,1:j-1));
    end
end
toc;
end

然后是调用函数的文件,C2P2.m:

clear;
clc;
in = -4:8/100:4;
tmp = 0.2:0.2:1.0;
fValue = [0.98 0.92 0.81 0.64 0.38];
tic;
out = NewtonPoly(tmp,fValue,in);
toc;
plot(in,out);

输出是:
Elapsed time is 0.000086 seconds.
Elapsed time is 0.001977 seconds.
Elapsed time is 0.002401 seconds.

其中第一个是第一个只有一层的循环所耗费的时间;第二个是两层循环所耗费的时间;最后的是总共的时间。第一个时间耗费的时间少是因为ptsLen比较短,所以不明显。

不过还是可以清晰地看出这个问题,for循环的效率太低。可以想象如果我在第一个循环中不矩阵化,时间肯定会耗费的更多一些。

矩阵化虽然比较好,不过想起来也比较费脑子,又要省空间又要省时间挺难办到的,感觉还是比较需要经验的。