一、horner求导数和积分
对于一个多项式P(x)，我们可以用horner方法（嵌套乘法）求出多项式的值，同样的，将其修改一下，我们也可以求出其积分和导数的值，积分时设常数项C=0：
function [ ans ] = hornerIx( arr, x )
%输入:
%行向量arr [an ...a0]
%变量x
%求多项式p(x)的积分值I(x)
count = length(arr);
ans = arr(1)/(count);
for i=(count-1):-1:1
ans = ans*x + arr(count - i+1)/i;
end
ans = ans*x;
end
function [ ans ] = hornerDx( arr, x )
%输入:
%行向量arr [an ...a0]
%变量x
%求多项式p(x)的导数值dp(x)
count = length(arr)-1;
ans = arr(1)*(count);
for i=count-1:-1:1
ans = ans*x + arr(count - i+1)*i;
end
end
二、拉格朗日逼近
1、poly()
poly()函数用于求以向量为解的方程或方阵的特征多项式。（第一个系数为1）
以方程为例，输入参数为一个行向量，表示多项式的根：[1,2,3,4]表示多项式的根为x=1、2、3、4，求解结果如下：
>> poly([1,2,3,4])
ans =
1 -10 35 -50 24
即多项式为$P(x) = x^4 -10x^3 + 35x^2 -50x +24$
2、conv()
conv（向量卷积运算）求两个向量卷积，简单理解其实就是多项式乘法。
进行多项式乘法时，输入参数为两个行向量conv([an,...a0],[bn...b0])
>> P = poly([1,-1])
P =
1 0 -1
>> Q = poly([2])
Q =
1 -2
>> conv(P,Q)
ans =
1 -2 -1 2
$conv(P, Q) = P(x)Q(x) = (x^2-1)(x-2) = x^3 -2x^2 - x +2$
3、拉格朗日逼近
给定N+1个点，得到N次拉格朗日多项式：
输入X、Y为长度为N+1的行向量，输出C为N次多项式系数，L为NxN的矩阵，矩阵第k行为$L_{Nk}()$的值。
function [ C, L ] = lagran( X,Y )
% input
% X、Y：N+1个点（x，y）
% output：
% C：
% L:
w = length(X);
n = w-1;
L = zeros(w,w);
for k = 1:n+1
V = 1;
for j = 1:n+1
if k~=j
V = conv(V,poly(X(j))) / (X(k)-X(j));
end
end
L(k,:) = V;
end
C = Y*L;
end
其中核心部分为V = conv(V,poly(X(j))) / (X(k)-X(j));， 这一个循环求的是拉格朗日系数多项式L，将每一个多项式存入矩阵L的每一行中，最后C = Y*L;求出拉格朗日多项式C。
下面是多个测例：
>> [C,L] = lagran([1,2,3],[1,2,3])
C =
0 1 0
L =
0.5000 -2.5000 3.0000
-1.0000 4.0000 -3.0000
0.5000 -1.5000 1.0000
>> [C,L] = lagran([-1,0],[-1,0])
C =
1 0
L =
-1 0
1 1
>> [C,L] = lagran([-1,0,1,2],[-1,0,1,8])
C =
1.0000 0 0.0000 0
L =
-0.1667 0.5000 -0.3333 0
0.5000 -1.0000 -0.5000 1.0000
-0.5000 0.5000 1.0000 0
0.1667 0 -0.1667 0
>> [C,L] = lagran([0,1,2],[0,1,8])
C =
3 -2 0
L =
0.5000 -1.5000 1.0000
-1.0000 2.0000 0
0.5000 -0.5000 0
三、牛顿插值多项式
给定N+1个点，得到N次牛顿插值多项式：
输入X、Y为长度为N+1的行向量，输出C为N次多项式系数，D为(N+1)x(N+1)的差商表。
function [ C, D ] = newpoly( X,Y )
% input
% X、Y：N+1个点（x，y）
% output：
% C：多项式系数
% D: 差商表
n = length(X);
L = zeros(n,n);
D(:,1) =Y'; % 将Y列向量存入D的第一列，初始化差商表f[xk]
% 求出差商表
for j = 2:n
for k = j:n
D(k,j) = (D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));
end
end
% 嵌套乘法求多项式C
C = D(n,n);
for k = (n-1):-1:1
C = conv(C,poly(X(k)));
m = length(C);
C(m) = C(m) + D(k,k);
end
end
>> [C,D] = newpoly([0,1,2,3,4],[0,0.75,2.25,3,2.25])
C =
0.0313 -0.4375 1.4688 -0.3125 0
D =
0 0 0 0 0
0.7500 0.7500 0 0 0
2.2500 1.5000 0.3750 0 0
3.0000 0.7500 -0.3750 -0.2500 0
2.2500 -0.7500 -0.7500 -0.1250 0.0313
>> [C,D] = newpoly([0,1,2,3],[0,0.75,2.25,3])
C =
-0.2500 1.1250 -0.1250 0
D =
0 0 0 0
0.7500 0.7500 0 0
2.2500 1.5000 0.3750 0
3.0000 0.7500 -0.3750 -0.2500
>> [C,D] = newpoly([4,5,6,7,8],[2, 2.23, 2.44, 2.64, 2.82])
C =
-0.0008 0.0200 -0.1842 0.9750 -0.0200
D =
2.0000 0 0 0 0
2.2300 0.2300 0 0 0
2.4400 0.2100 -0.0100 0 0
2.6400 0.2000 -0.0050 0.0017 0
2.8200 0.1800 -0.0100 -0.0017 -0.0008
