2014年8月21日 星期四

數值積分

許多函數的積分結果並沒有解析解,此時就要利用數值積分算出近似值,常用的方法有
  • 矩形法(維基百科 中文 en
  • 梯形法(維基百科 中文 en
  • 辛普森1/3法則及3/8法則(維基百科 中文 en

為了完成這些計算,我們可以自己寫程式,也可以利用LibreOffice Calc、Microsoft Excel這類的工作表軟體,甚至還可以利用Wolfram Alpha網站(http://www.wolframalpha.com/)查詢答案。以運算速度以及可以分割的區間數量而言,自己寫程式一定是最好的。

以下是我用LibreOffice Calc製作的檔案:數值積分.ods
數值積分_calc_pic4.png


另外我也用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即可。但如果想要看到進一步的資料,就必須付費註冊才行。




2 則留言:

  1. 回覆
    1. 沒想到這麼久以前的文章會被找到,希望這篇文章能幫助到你。

      刪除