牽著老婆滿街逛

嚴以律己,寬以待人. 三思而后行.
GMail/GTalk: yanglinbo#google.com;
MSN/Email: tx7do#yahoo.com.cn;
QQ: 3 0 3 3 9 6 9 2 0 .

Blitz++ 計算二階導數

// 整理?by?RobinKin

#include? < blitz / vector.h >

using ? namespace ?blitz;

int ?main()

{

???? // ?In?this?example,?the?function?cos(x)^2?and?its?second?derivative

???? // ?2?(sin(x)^2?-?cos(x)^2)?are?sampled?over?the?range?[0,1).

???? // ?The?second?derivative?is?approximated?numerically?using?a

???? // ?[?1?-2??1?]?mask,?and?the?approximation?error?is?computed.

/* ?cos(x)^2??的二階導數?是2?(sin(x)^2?-?cos(x)^2)

下面在[0,1)?的范圍中，以delta為步長?，用[?1?-2??1?]?的方法計算二階導數

看看和精確值的誤差

???? const ? int ?numSamples? = ? 100 ;?????????????? // ?Number?of?samples

???? double ?delta? = ? 1 .? / ?numSamples;?????????? // ?Spacing?of?samples

????Range?R( 0 ,?numSamples? - ? 1 );?????????????? // ?Index?set?of?the?vector

???? // ?Sample?the?function?y?=?cos(x)^2?over?[0,1)

???? //

???? // ?An?object?of?type?Range?can?be?treated?as?a?vector,?and?used

???? // ?as?a?term?in?vector?expressions.

???? //

???? // ?The?initialization?for?y?(below)?will?be?translated?via?expression

???? // ?templates?into?something?of?the?flavour

???? //

???? // ?for?(unsigned?i=0;?i?<?99;?++i)

???? // ?{

???? // ?????double?_t1?=?cos(i?*?delta);

???? // ?????y[i]?=?_t1?*?_t1;

???? // ?}

????

????Vector < double > ?y? = ?sqr(cos(R? * ?delta));

???? // ?Sample?the?exact?second?derivative

????Vector < double > ?y2exact? = ? 2.0 ? * ?(sqr(sin(R? * ?delta))? - ?sqr(cos(R? * ?delta)));

???? // ?Approximate?the?2nd?derivative?using?a?[?1?-2??1?]?mask

???? // ?We?can?only?apply?this?mask?to?the?elements?1?..?98,?since

???? // ?we?need?one?element?on?either?side?to?apply?the?mask.

????Range?I( 1 ,numSamples - 2 );

????Vector < double > ?y2(numSamples);

????y2(I)? = ?(y(I - 1 )? - ? 2 ? * ?y(I)? + ?y(I + 1 ))? / ?(delta * delta);

???? // ?The?above?difference?equation?will?be?transformed?into

???? // ?something?along?the?lines?of

???? //

???? // ?double?_t2?=?delta*delta;

???? // ?for?(int?i=1;?i?<?99;?++i)

???? // ?????y2[i]?=?(y[i-1]?-?2?*?y[i]?+?y[i+1])?/?_t2;

???? // ?Now?calculate?the?root?mean?square?approximation?error:

???? double ?error? = ?sqrt(mean(sqr(y2(I)? - ?y2exact(I))));

???? // ?Display?a?few?elements?from?the?vectors.

???? // ?This?range?constructor?means?elements?1?to?91?in?increments

???? // ?of?15.

????Range?displayRange( 1 ,? 91 ,? 15 );

????cout? << ? " Exact?derivative: " ? << ?y2exact(displayRange)? << ?endl

????????? << ? " Approximation:??? " ? << ?y2(Range(displayRange))? << ?endl

????????? << ? " RMS?Error:??????? " ? << ?error? << ?endl;

???? return ? 0 ;

}

Output:

Exact?derivative:[???? - 1.9996 ?? - 1.89847 ?? - 1.62776 ?? - 1.21164 ? - 0.687291 ? - 0.1015495

???]

Approximation:???[??? - 1.99953 ?? - 1.89841 ??? - 1.6277 ??? - 1.2116 ? - 0.687269 ? - 0.1015468

???]

RMS?Error:??????? 4.24826e-05

posted on 2006-07-01 13:23 楊粼波閱讀(699) 評論(0) 編輯收藏引用

只有注冊用戶登錄后才能發表評論。
【推薦】100%開源！大型工業跨平臺軟件C++源碼提供，建模，組態！



網站導航: 博客園 IT新聞 BlogJava 博問 Chat2DB 管理

牽著老婆滿街逛

導航

統計

公告

常用鏈接

留言簿(11)

隨筆分類(466)

隨筆檔案(1513)

文章分類(46)

文章檔案(45)

相冊

收藏夾(39)

工具官網

技術網站

開源網站

其他窩點

收藏網站

銀行官網

友情鏈接

資源共享

搜索

積分與排名

最新評論

閱讀排行榜

Blitz++ 計算二階導數