多项式插值算法 python实现

算法思想：用多项式逼近原函数

import numpy as np
from numpy import *
import numpy.linalg as lg    #numpy的线性代数函数库 linalg
import math
x=[]
y=[]
N=10
pi=math.pi
#形成十个(x,y）点
for i in range(N):
    x.append(round((-1+(2/N)*(N-i)),3))
def function(x1):
    return math.sin(pi*x1)
for i in range(len(x)):
    y.append(function(x[i]))


n = len(x)
X= zeros((n, n))
for i in range(n):
    for j in range(n):
        X[i][j] = math.pow(x[i],j)
print(X)    #范德蒙
A=[]    #系数矩阵
XT=lg.inv(X)  #求X得逆矩阵
Y=zeros((n,1))
for i in range(len(y)):
    Y[i][0]=y[i] 

A=lg.solve(X,Y)    #解系数矩阵
def polynomial_interploate(x1):
    P=0
    for i in range(len(x)):
        P+=A[i][0]*math.pow(x1,i)
    return P

function y = lagrange(x0, y0, x)
    %x0，y0是插值节点
    %x是待计算其它点，y为这些点上的函数值
    n=length(x0);
    m=length(x);
    for i = 1:m
        z = x(i);
        y(i) = baseLagrange(x0, y0, z, n);
    end
%嵌套插值基函数
function s = baseLagrange(x0, y0, z, n)
    s = 0.0;
    for k = 1:n
        p = 1.0;
        for j = 1:n
            if j ~= k
                p = p*(z - x0(j))/(x0(k) - x0(j));
            end
        end
        s = p*y0(k) + s;
    end
end
end

#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double type;
type lagrange(vector<type> x,vector<type> y,type X)
{
	int n=x.size()-1;
	type ans=0;
	rep(k,0,n)
	{
		type temp=y[k];
		rep(j,0,n)if(j!=k)temp*=(X-x[j])/(x[k]-x[j]);
		ans+=temp;
	}
	return ans;
}
int main()
{
	vector<type> x;
	vector<type> y;
	rep(i,0,4)
        x.push_back(i),
        y.push_back(1.0*i*i+1.0*i/6);
	type X;
	while(cin>>X)cout<<lagrange(x,y,X)<<endl;
	return 0;
}
/*
f(x) = x*x + x/6
    0.5 -> 0.333333
    1.2 -> 1.64
*/

#define rep(i,a,b) for(LL i=a;i<=b;i++)
using namespace std;
typedef long long LL;
typedef long long type;
const LL mod=1e9+7;
LL Pow(LL a,LL b,LL mod){
    LL res=1;
    while(b>0){
        if(b&1)res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}
LL inv[10][10];
type lagrange(vector<type> x,vector<type> y,type X)
{
	int n=x.size()-1;
	type ans=0;
	rep(i,0,n)
	{
		type temp=y[i];
		rep(j,0,n)if(j!=i)temp=(X-x[j])*inv[i][j]%mod*temp%mod;
		ans+=temp;
	}
	return (ans%mod+mod)%mod;
}
int main()
{
	vector<type> x={0,1,2,3,4};
	vector<type> y={0,1,5,15,35};
	int n=x.size()-1;
	rep(i,0,n)
        rep(j,0,n)
            inv[i][j]=Pow(x[i]-x[j],mod-2,mod );
	type X;
	int t;scanf("%d",&t);
	while(t--)scanf("%lld",&X),printf("%lld\n",lagrange(x,y,X));
	return 0;
}