随笔 - 51, 文章 - 1, 评论 - 41, 引用 - 0
数据加载中……

求向量的欧拉范数

      求欧拉范数是一个比较简单的算法，似乎没有什么可说的，一般代码如下：
/** \fn double enorm1(long n, const double* x)
* \brief 求欧拉范数，简单算法。
* \param [in] n 向量长度
* \param [in] x 向量值 
* \return 欧拉范数
*/
double enorm1(long n, const double* x)
{
    double ret = 0.0;
    long i = 0;
    for (i=0; i<n; ++i)
        ret += x[i]*x[i];
        
    return sqrt(ret);
}

然而近日在学习minpack时，发现它求欧拉范数的函数enorm则要复杂得多。仔细比较发现上面算法有不足，它没有考虑溢出的情况，当x[i]为很小或者很大的数时，x[i]*x[i]则会下溢或者上溢，最后结果可能不准确。但x[i]*x[i]只是中间结果，最后的范数数量级应该和x[i]相同。它溢出的可能性要小得多。改进中间结果则可以改进算法。

/** \fn double enorm2(long n, const double* x)
* \brief 求欧拉范数，考虑溢出的算法。
* \param [in] n 向量长度
* \param [in] x 向量值 
* \return 欧拉范数
*/
double enorm2(long n, const double *x)
{
    double ret = 0.0;
    long i= 0;
    double xmax = 0.0;
    
    for (i=0; i<n; ++i)
    {
        double xabs = fabs(x[i]);
        
        /* 这个比较方式不需要考虑xabs和xmax为0的情况 */
        if (xabs < xmax)
        {
            ret += (xabs/xmax)*(xabs/xmax);
        }
        else if (xabs == xmax)  
        {
            ret += 1;
        }
        else
        {
            ret = 1 + ret*(xmax/xabs)*(xmax/xabs);
            xmax = xabs;
        }
    }
    return sqrt(ret) * xmax;
}

   将中间结果改成(x[i]/xmax)*(x[i]/xmax)，降低了溢出的可能。当然该算法没有区分可能溢出和不溢出的数，计算量较大。下面的算法是仿照minpack的enorm函数编写：

/** \fn double enorm3(long n, const double* x)
* \brief 求欧拉范数，仿照minpack的enorm函数
* \param [in] n 向量长度
* \param [in] x 向量值 
* \return 欧拉范数
*/
double enorm3(long n, const double *x)
{
    double s1 = 0.0;    
    double s2 = 0.0;    
    double s3 = 0.0;    
    
    /* 上溢和下溢的边界，不一定要十分精确 */
    const double dwarf = 1.483e-154;    /* 下溢的边界 */
    const double giant = 1.341e154 / n; /* 上溢的边界 */
    
    double x1max = 0.0; 
    double x3max = 0.0; 
    
    long i = 0;
    for (i=0; i<n; ++i)
    {
        double xabs = fabs(x[i]);
        
        if (xabs > dwarf && xabs < giant)
        {
            s2 += xabs*xabs;
        }
        else if (xabs <= dwarf)
        {
            /* 这个比较方式不需要考虑xabs和xmax为0的情况 */
            if (xabs < x3max)
            {
                s3 += (xabs/x3max)*(xabs/x3max);
            }
            else if (xabs == x3max)
            {
                s3 += 1;
            }
            else
            {
                s3 =1 + s3*(x3max/xabs)*(x3max/xabs);
                x3max = xabs;
            }
        }
        else /* if (xabs >= giant) */
        {
            /* 不需要考虑xabs和xmax为0的情况 */
            if (xabs <= x1max)
            {
                s1 += (xabs/x1max)*(xabs/x1max);
            }
            else if (xabs == x1max)
            {
                s1 += 1;
            }
            else
            {
                s1 =1 + s1*(x1max/xabs)*(x1max/xabs);
                x1max = xabs;
            }
        }
    }
    
    if (s1 != 0.0)
    {
        return x1max*sqrt(s1 + s2/x1max/x1max);
    }
    else if (s2 != 0.0)
    {
        return sqrt(s2 + x3max*x3max*s3);
        /* 下面为minpack中enorm的代码，好像没有必要 
        if (s2 >= x3max)
            return sqrt(s2 * (1 + x3max/s2*x3max*s3));
        else
            return sqrt(x3max * (s2/x3max + x3max*s2));
        */
    }
    else
    {
        return x3max * sqrt(s3);
    }
}

 

下面是测试结果：

<td style="BORDER-RIGHT: rgb(0,0,0) 0.5pt solid; PADDING-RIGHT: 5.4pt; BORDER-TOP: medium none; PADDING-LEFT: 5.4pt; PADDING-BOTTOM: 0pt; BO

 

posted on 2009-05-07 13:07 lemene 阅读(1508) 评论(0)  编辑 收藏 引用

向量

enorm1

enorm2

enorm3

{1, 2, 3}

3.74166

3.74166

3.74166

{1e200, 2e200, 3e200}

1.#INF