http://keejko.blogspot.com/2014/08/blog-post.html
數值積分
許多函數的積分結果並沒有解析解，此時就要利用數值積分算出近似值，常用的方法有
為了完成這些計算，我們可以自己寫程式，也可以利用LibreOffice Calc、Microsoft Excel這類的工作表軟體，甚至還可以利用Wolfram Alpha網站(http://www.wolframalpha.com/)查詢答案。以運算速度以及可以分割的區間數量而言，自己寫程式一定是最好的。
以下是我用LibreOffice Calc製作的檔案：數值積分.ods
另外我也用c語言寫了幾個小程式，主要是從維基百科上找到的程式碼改寫的，測試的環境為Windows 8.1 64位元版本、Dev-C++ 5.7.1版，以下是原始碼：
/* 數值積分：梯形法
由維基百科上的程式碼改寫而成
來源： http://zh.wikipedia.org/wiki/矩形法
改作者：王一哲 Yi-Zhe Wang
日期：Aug. 21, 2014 */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// 定義函數f(x)，也可以回傳其它數學函數
double f(double x){
return sin(x);
}
// 定義運算式
double trapezoidal_integrate(double a, double b, int subintervals){
double result;
double interval;
int i;
interval=(b-a)/subintervals;
result=0;
//首項、末項係數為1，其餘係數為2
for(i=1;i<=subintervals;i++){
if((i==1)||(i==subintervals)){
result+=f(a+interval*(i-1));
} else {
result+=2*f(a+interval*(i-1));
}
}
result*=interval/2;
return result;
}
//主程式
int main(void){
double integral;
integral=trapezoidal_integrate(0,2,500000);
printf("Integral: %.15f \n",integral);
system("pause");
return 0;
}
/* 數值積分：辛普森1/3法則
由維基百科上的程式碼改寫而成
來源： http://zh.wikipedia.org/wiki/矩形法
改作者：王一哲 Yi-Zhe Wang
日期：Aug. 21, 2014 */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// 定義函數f(x)，也可以回傳其它數學函數
double f(double x){
return sin(x);
}
// 定義運算式
double simpson1_3_integrate(double a, double b, int subintervals){
double result;
double interval;
int i;
interval=(b-a)/subintervals;
result=0;
//首項、末項係數為1，其餘每2項1組，第1項係數為4，第2項係數為2
for(i=0;i<=subintervals;i++){
if((i==0)||(i==subintervals)){
result+=f(a+interval*i);
} else if((i%2)==1){
result+=4*f(a+interval*i);
} else if((i%2)==0){
result+=2*f(a+interval*i);
}
}
result*=interval/3;
return result;
}
//主程式，分割的區間數只能是偶數
int main(void){
double integral;
integral=simpson1_3_integrate(0,2,500000);
printf("Integral: %.15f \n",integral);
system("pause");
return 0;
}
/* 數值積分：辛普森3/8法則
由維基百科上的程式碼改寫而成
來源： http://zh.wikipedia.org/wiki/矩形法
改作者：王一哲 Yi-Zhe Wang
日期：Aug. 21, 2014 */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// 定義函數f(x)，也可以回傳其它數學函數
double f(double x){
return sin(x);
}
// 定義運算式
double simpson3_8_integrate(double a, double b, int subintervals){
double result;
double interval;
int i;
interval=(b-a)/subintervals;
result=0;
//首項、末項係數為1，其餘每3項1組，1、2項係數為3，第3項係數為2
for(i=0;i<=subintervals;i++){
if((i==0)||(i==subintervals)){
result+=f(a+interval*i);
} else if(((i%3)==1)||((i%3)==2)){
result+=3*f(a+interval*i);
} else if((i%3)==0){
result+=2*f(a+interval*i);
}
}
result*=interval*3/8;
return result;
}
//主程式，分割的區間數只能是3的倍數
int main(void){
double integral;
integral=simpson3_8_integrate(0,2,500001);
printf("Integral: %.15f \n",integral);
system("pause");
return 0;
}
如果要利用Wolfram Alpha網站查詢答案，就在方框內輸入
Integrate[f(x),{x,a,b}]
上式中的f(x)請換成要積分的函數，a換成積分下限，b換成積分下限，最後按下Enter即可。但如果想要看到進一步的資料，就必須付費註冊才行。
沒有留言:
張貼留言