求歐拉范數是一個比較簡單的算法，似乎沒有什么可說的，一般代碼如下：
/** \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);
}
}
下面是測試結果：
|
向量
|
enorm1
|
enorm2
|
enorm3
|
|
{1, 2, 3}
|
3.74166
|
3.74166
|
3.74166
|
|
{1e200, 2e200, 3e200}
|
1.#INF
|
<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
導航
|
|
| 26 | 27 | 28 | 29 | 30 | 1 | 2 |
| 3 | 4 | 5 | 6 | 7 | 8 | 9 |
| 10 | 11 | 12 | 13 | 14 | 15 | 16 |
| 17 | 18 | 19 | 20 | 21 | 22 | 23 |
| 24 | 25 | 26 | 27 | 28 | 29 | 30 |
| 31 | 1 | 2 | 3 | 4 | 5 | 6 |
常用鏈接
留言簿(4)
隨筆檔案
文章檔案
搜索
最新隨筆
最新評論
閱讀排行榜
評論排行榜
青青草原综合久久大伊人导航_色综合久久天天综合_日日噜噜夜夜狠狠久久丁香五月_热久久这里只有精品