[問題] numerical recipe C++ 積分的問題

看板C_and_CPP (C/C++)作者peter74123時間16年前 (2010/01/21 14:49)推噓0(0推 0噓 1→)

留言1則, 1人參與討論串1/1

請求大家幫看一下這個問題拜託了... 我是理工科系碩二生但是程式只會最基本的 numerical recipe我只會"使用" 但是運作原裡完全不懂... 我碰到的問題是這樣的我使用numerial recipe第二版的qsimp這個routine來做數值積分所以我可以照著numerical recipe書上的作法定義一個自訂函數func 此函數就是f(x) 然後在main裡面NR::qsimp(func,下界,上界); 算出f(x)的積分但是我希望此函數可以是f(x;a,b,c,d...)的形式也就是說改成可以輸入a b c d等等參數的數值後做積分如果我只算一次當然就沒差但是實際上我是要算幾百次每次的a b c d數值都不同我有試過改函數的形式例如DP func(const DP x) 改成 DP func(const DP x, DP a, DP b) 或是用其他方法去改函式但是就是不行可能因為func,qsimp等等都是定義在nr內的東西吧所以要到nr裡面改嘛？應該有其他辦法吧.... 希望有使用過numerical recipe的版友可以幫我解決這問題感激不進！程式碼 --- #include <iostream> #include <iomanip> #include <cmath> #include <cstdlib> #include "NR/nr.h" using namespace std; // Driver for routine trapzd DP NR::trapzd(DP func(const DP), const DP a, const DP b, const int n) { DP x,tnm,sum,del; static DP s; int it,j; if (n==1){ return (s=0.5*(b-a)*(func(a)+func(b))); } else { for (it=1,j=1;j<n-1;j++) it <<= 1; tnm=it; del=(b-a)/tnm; x=a+0.5*del; for (sum=0.0,j=0;j<it;j++,x+=del) sum += func(x); s=0.5*(s+(b-a)*sum/tnm); return s; } } // Driver for routine qsimp DP NR::qsimp(DP func(const DP), const DP a, const DP b) { const int JMAX=20; const DP EPS=1.0e-10; int j; DP s,st,ost=0.0,os=0.0; for (j=0;j<JMAX;j++) { st=trapzd(func,a,b,j+1); s=(4.0*st-ost)/3.0; if (j>5) if (fabs(s-os) < EPS*fabs(os) || (s==0.0 && os==0.0)) return s; os=s; ost=st; } nrerror("Too many steps inroutine qtrap"); return 0.0; } // 自訂函數 DP func(const DP x) { //想把這四個參數改成可調整，例如用輸入的 DP rS = 0.194951; DP rhoS = 9.2121e+15; DP R = 0.8; DP rsub = 0.2; DP critical= 3.524689e+15; DP r = sqrt(R*R+rsub*rsub+2.0*R*rsub*cos(x)); return (1.0/critical)*2.0*rhoS*rS*(1.0/(1.0-r*r/(rS*rS)))*(-1.0+2.0*atan(sqrt((r/rS-1.0)/(r/rS+1.0)))/sqrt(r*r/(rS*rS)-1.0)); } int main(void) { double pi=4.0*atan(1.0); double s; s=NR::qsimp(func,0,2.0*pi); s=s/(2.0*pi); cout << fixed << setprecision(6); cout<<s<<endl; return 0; } -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 140.112.4.187 ※ 編輯: peter74123 來自: 140.112.4.187 (01/21 14:56)