例題4-2梯行法(Trapezoidal method)計算 exp( - (z^2)/2 ) / (2* Pi ) 在[-3.5, -3.0]的定積分

<C++程式語言>例題4-2計算 exp( - (z^2)/2 ) / (2* Pi ) 在[-3.5, -3.0]的定積分

// C++ program to implement Trapezoidal rule
#include<stdio.h>
#include<math.h>
# define M_PI           3.14159265358979323846  /* pi */
 
// A sample function whose definite integral's
// approximate value is computed using Trapezoidal
// rule
float y(float x)
{
    // Declaring the function f(x) = 1/(1+x*x)
    return (1/sqrt(2*M_PI))*exp(-pow(x,2)/2);
}

float f(float x)
{
    // Declaring the function f(x) = 1/(1+x*x)
    return (pow(x,2)-1)*exp(-pow(x,2)/2);
}

 
// Function to evalute the value of integral
float trapezoidal(float a, float b, float n)
{
    // Grid spacing
    float h = (b-a)/n;
 
    // Computing sum of first and last terms
    // in above formula
    float s = y(a)+y(b);
 
    // Adding middle terms in above formula
    for (int i = 1; i < n; i++)
        s += 2*y(a+i*h);
 
    // h/2 indicates (b-a)/2n. Multiplying h/2
    // with s.
    return (h/2)*s;
}
 
// Driver program to test above function
int main()
{
    // Range of definite integral
    float x0 = -3.5;
    float xn = -3.0;
 
    // Number of grids. Higher value means
    // more accuracy
    int n = 10;
    float ans1=trapezoidal(x0, xn, n);
    //float ans2=f(xn)-f(x0);
    //printf("Value of integral is %6.7f\n",f(xn));
    //printf("Value of integral is %6.7f\n",f(x0));
   
    printf("Value of integral is %6.7f\n",ans1);
       
   
    return 0;
}

輸出畫面
Value of integral is 0.0011194