注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

民主与科学

独立之人格,自由之思想

 
 
 

日志

 
 

拉格朗日插值法  

2010-09-14 15:19:45|  分类: 数学相关 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
拉格朗日插值法
拉格朗日插值法可以帮助我们解决以下的问题
已知x取值0,1,-1,2时, f{x}取值2,2,0,6 
求x=3时f{x}的值。

示例1
        int xs[]={0,1,-1,2};
        int ys[]={2,2,0,6};
        
//f(3)?
        int val=lagrangePolynomial(3,xs,ys);
        
    static int  lagrangePolynomial(int x, int xs[],int ys[]) 
    {
        int c1,c2; 
        int i,j,n; 
        n=xs.length;
        c1=0; 
        for(i=0;i<n;i++){ 
            c2=ys[i]
            for(j=0;j<n;j++)
            {
                if(j!=i)
                    c2=c2*(x-xs[j])/(xs-xs[j]); 
            }
            c1+=c2; 
        } 
       
 return c1;
    }

Lagrange插值方法的核心就是构造一组基函数
如果插值点是{(x_i,y_i)}i=1..n,那么希望构造出一组多项式F_i(x)使得
F_i(x_i)=1, F_i(x_j)=0 (j!=i)
其实就是构造一个函数, 这个函数在其中一点的值为1, 其它点的值为0 。 
这样的话把n个这样的函数加权加起来得到的函数就是在每个点上的值都是需要的了。
注意1这里的多项式是一组。函数即是指多项式。其实针对每个插入点构造一个多项式(函数)
注意2:关于如何得到多项式的请参阅http://hwiechern.blog.163.com/blog/static/106796622007913104859712/
也就是说要构造“只受其中一个点影响”(这种讲法比较粗糙,因为和其他点的位置还是有关系)的函数。
如果这一点能办到,那么只要取f(x)=sum(y_i*F_i(x))就是所要的插值多项式。
Lagrange的插值方法其实就是直接构造出上述基函数:
F_i(x) = prod(x-x_j) / prod(x_i-x_j),其中prod是关于所有不等于i的j求乘积,
直接就可以验证F_i(x)满足前面提到的条件,
因为分子相当于确定了F_i(x)的所有根,分母则是归一化系数。


对于一些很复杂的函数运算。也可以用拉格朗日插值法来计算以提高效率
我们可以先用函数运算均匀的生成一些x点和其对应得数值对。
然后就可以用拉格朗日插值法来计算以提高效率。
具体使用见实例1
    static void lagrangeDemo()
    {
        int k=10000;
        int a=2;
        int x=0;
        int y=0;
        int m=0;
        float data[]=new float[k];
        Random random=new Random(System.currentTimeMillis());
        for(int i=0;i<k;i++)
        {
            do{
            m=random.nextInt((int)(kMaxX-kMinX));
            data=m+1/random.nextInt();
            }while(data<kMinX||data>kMaxX);
        }
        float res0,res1;
        float diff=0;
        float diffMax=0;
        float diffRate=0;
        float diffRateMax=0;
        float diffAll=0;
        for(int i=0;i<k;i++)
        {
            res0=count(a,data);
            res1=function(a,data);
            diff=(res1-res0);
            if(Math.abs(diff)>diffMax)
                diffMax=diff;
            diffAll+=diff;
            if(res1!=0)
            {
                diffRate=diff*100/res1;
                if(Math.abs(diffRate)>diffRateMax)
                    diffRateMax=diffRate;
            }
        }
        System.out.println("diffAll:"+diffAll+"diffMax:"+(diffMax)+"diffRateMax:"+diffRateMax+"%");
        long time=System.currentTimeMillis();
        for(int i=0;i<k;i++)
        {
            res0=count(a,data);
        }
        System.out.println(" lagrange time:"+(System.currentTimeMillis()-time));
        time=System.currentTimeMillis();
        for(int i=0;i<k;i++)
        {
            res1=function(a,data);
        }
        System.out.println("Normal time:"+(System.currentTimeMillis()-time));
    }
    static HashMap <String,float[][]>hashMap=new HashMap<String,float[][]>();
    final static int kSection=100;
    static float count(int a,float x)
    {
        String key=""+a;
        float f[][]=hashMap.get(key);
       
if(f==null)
       {
                f=getPairs(a);
                hashMap.put(key, f);
        }
        
return lagrangePolynomial(x,f[0],f[1]);
    }
    final static float kMinX=0;
    final static float kMaxX=100;
   
 static float[][] getPairs(int a)
    {
        float pairs[][]=new float[2][kSection+1];
        float xs[]=pairs[0];
        float ys[]=pairs[1];
        int n=kSection;
        
for(int i=0;i<=n;i++)
        {
            xs=kMinX+(kMaxX-kMinX)*i/n;
            ys=function(a,xs);
        }
        
return pairs;
    }
    
static float function(int a,float x)
    {
        float y=x*x*x*x/(a*a+x*x*x*x+x*x);
        int n=1000;
        int i,j;
        
for(i=0;i<n;i++)
        {
            
for(j=0;j<n;j++)
            {
                y=y*n*(n+2)/((n+1)*(n+3))-1;
            }
        }
        return y;
    }
    
static int  lagrangePolynomial(int x, int xs[],int ys[]) 
    {
        int c1,c2; 
        int i,j,n; 
        n=xs.length;
        c1=0; 
        
for(i=0;i<n;i++){ 
            c2=ys
            
for(j=0;j<n;j++)
            {
                if(j!=i)
                    c2=c2*(x-xs[j])/(xs-xs[j]); 
            }
            c1+=c2; 
        } 
        
return c1;
    }
    
static float  lagrangePolynomial(float x, float xs[],float ys[]) 
    {
        float c1,c2; 
        int i,j,n; 
        n=xs.length;
        c1=0; 
        
for(i=0;i<n;i++){ 
            c2=ys
            
for(j=0;j<n;j++)
            {
                if(j!=i)
                    c2=c2*(x-xs[j])/(xs-xs[j]); 
            }
            c1+=c2; 
        } 
       
 return c1;
    }

运行结果
diffAll:NaNdiffMax:4.8828125E-4diffRateMax:-6.088556E-6%
lagrange time:687
Normal time:99359

注意:只有函数f(x)非常复杂,且f(x)的值随X变化不是很大时,用拉格朗日插值法来计算才有意义。
否则拉格朗日插值法的效率比直接用函数计算还要低。
关于拉格朗日插值法的更多数学知识请参考
http://hwiechern.blog.163.com/blog/static/106796622007913104859712/
  评论这张
 
阅读(613)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017