一個地球人在臺北: 數值積分

許多函數的積分結果並沒有解析解，此時就要利用數值積分算出近似值，常用的方法有

矩形法（維基百科　中文　en）
梯形法（維基百科　中文　en）
辛普森1/3法則及3/8法則（維基百科　中文　en）

為了完成這些計算，我們可以自己寫程式，也可以利用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 <stdlib.h>
#include <math.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 <stdlib.h>
#include <math.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 <stdlib.h>
#include <math.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即可。但如果想要看到進一步的資料，就必須付費註冊才行。