都是比較復雜的函數,可以先借助help、type和edit命令查看函數源代碼;數學三大庫GNU的GSL,blitz++以及MTL對積分也有支持;
要實現像matlab壹樣輸入靈活的函數,個人感覺還是很不容易的,但是面向具體的業務問題,在函數及條件給定前提下獲得符號積分、數值積分近似解的實現還不是很難的,因為流程不外乎就是矩形法、梯形法和拋物線法,eg:用梯形法(1000等份)求以下三個函數的定積分(當然如果作為工業使用這個for循環需要解決累加溢出和性能問題):
f1(x)=1+x*x(積分區間0到1)
f2(x)=1+x+x*x+x*x*x(積分區間0到2)
f3(x)=x/(1+x*x)(積分區間0到3.5)
#include?"iostream.h"#define?N?1000
double?f1(double?x)//f1(x)=1+x*x
{
return(1+x*x);}
double?f2(double?x)//f2(x)=1+x+x*x+x*x*x
{
return(1+x+x*x+x*x*x);
}
double?f3(double?x)//f3(x)=x/(1+x*x)
{
return(x/(1+x*x));
}
double?integral(double(*?fun)(double),double?a,double?b)//求fun在a,b區間定積分
{
double?s,h;
int?i;
s=((*?fun)(a)+(*?fun)(b))/2.0;
h=(b-a)/N;
for(i=1;i<N;i++)
s+=(*fun)(a+i*h);return(s*h);
}
main()
{
double?y1,y2,y3;
y1=integral(f1,0.0,1.0);
y2=integral(f2,0.0,2.0);
y3=integral(f3,0.0,3.5);
cout<<"y1=?"<<y1<<";?y2=?"<<y2<<";?y3=?"<<y3<<endl;
return?0;
}還可以參考課本上的這些積分算法:
常用數值積分算法實現:
1 梯形求積法(TRAPZD->QTRAP) 2 辛普森(Simpson)求積法(QSIMP) 3 龍貝格(Romberg)求積法(QROMB) 4 反常積分(MIDPNT, MIDINF, MIDSQL, MIDEXP->QROMO) 5 高斯求積法(QGAUS, GAULEG) 6 三重積分(QUAD3D)