習題 EX4-8-2 梯行法雙重積分---C語言

//習題 EX4-8-2 梯行法雙重積分

/* 習題 8-2.c based on Trapezoidal Rule to
 * compute the double integral.
 */
#include <stdio.h>
#include <math.h>
#define F(x,y)  ( sin (x+y))
#define C(x)    (log(x))
#define D(x)    (3+ exp(0.2*x))
void gy(int);
double x[50],y[50][50],gg[50],hy[50];
void main()
{
   int i,j,n;
   double a,b,hx=0.0,sum=0.0,ts=0.0;
   n=10;
   a=1.0;
   b=3.0;
 
   hx=(b-a)/n;
   for(i=0;i<=n;i++)
   {
      x[i]=a+i*hx;
      hy[i]=(D(x[i])-C(x[i]))/n;
      for(j=0;j<=n;j++)
         y[i][j]=C(x[i])+j*hy[i];
   }
   gy(n);
   for(i=1;i<=n-1;i++)
      sum=sum+gg[i];
   ts=(hx/2.0)*(gg[0]+gg[n]+2*sum);
   printf("T%d=%.6lf\n",n,ts);
   return;
}

void gy(int n)
{
   int i,j;
   double sum=0.0;
   for(i=0;i<=n;i++)
   {
      for(j=1;j<=n-1;j++)
      {
  sum=sum+F(x[i],y[i][j]);
      }
      gg[i]=(hy[i]/2.0)*(F(x[i],y[i][0])+
     F(x[i],y[i][n])+2*sum);
      sum=0.0;
   }
   return;
}

輸出結果
T10=-2.562917

習題 EX4-8-3 梯行法雙重積分 -----C語言

習題 EX4-8-3 梯行法雙重積分 C語言

/* 習題 8-3.c based on Trapezoidal Rule to
 * compute the double integral.
 */
#include <stdio.h>
#include <math.h>
#define F(x,y)  (4*x*y)
#define C(x)    (0)
#define D(x)    (1- pow(x,2))
void gy(int);
double x[50],y[50][50],gg[50],hy[50];
void main()
{
   int i,j,n;
   double a,b,hx=0.0,sum=0.0,ts=0.0;
   n=10;
   a=0.0;
   b=1.0;
 
   hx=(b-a)/n;
   for(i=0;i<=n;i++)
   {
      x[i]=a+i*hx;
      hy[i]=(D(x[i])-C(x[i]))/n;
      for(j=0;j<=n;j++)
         y[i][j]=C(x[i])+j*hy[i];
   }
   gy(n);
   for(i=1;i<=n-1;i++)
      sum=sum+gg[i];
   ts=(hx/2.0)*(gg[0]+gg[n]+2*sum);
   printf("T%d=%.6lf\n",n,ts);
   return;
}

void gy(int n)
{
   int i,j;
   double sum=0.0;
   for(i=0;i<=n;i++)
   {
      for(j=1;j<=n-1;j++)
      {
  sum=sum+F(x[i],y[i][j]);
      }
      gg[i]=(hy[i]/2.0)*(F(x[i],y[i][0])+
     F(x[i],y[i][n])+2*sum);
      sum=0.0;
   }
   return;
}

輸出結果
T10=0.331650

習題 EX4-8-1 梯行法雙重積分----C語言

習題 EX4-8-1 梯行法雙重積分----C語言

/* 習題 8-1.c based on Trapezoidal Rule to
 * compute the double integral.
 */
#include <stdio.h>
#include <math.h>
#define F(x,y)  (x+y)
#define C(x)    (pow(x,2))
#define D(x)    (x+2)
void gy(int);
double x[50],y[50][50],gg[50],hy[50];
void main()
{
   int i,j,n;
   double a,b,hx,sum=0.0,ts;
   n=10;
   a=-1.0;
   b=2.0;
 
   hx=(b-a)/n;
   for(i=0;i<=n;i++)
   {
      x[i]=a+i*hx;
      hy[i]=(D(x[i])-C(x[i]))/n;
      for(j=0;j<=n;j++)
         y[i][j]=C(x[i])+j*hy[i];
   }
   gy(n);
   for(i=1;i<=n-1;i++)
      sum=sum+gg[i];
   ts=(hx/2.0)*(gg[0]+gg[n]+2*sum);
   printf("T%d=%.6lf\n",n,ts);
   return;
}
void gy(int n)
{
   int i,j;
   double sum=0.0;
   for(i=0;i<=n;i++)
   {
      for(j=1;j<=n-1;j++)
      {
  sum=sum+F(x[i],y[i][j]);
      }
      gg[i]=(hy[i]/2.0)*(F(x[i],y[i][0])+
     F(x[i],y[i][n])+2*sum);
      sum=0.0;
   }
   return;
}


輸出結果
T10=9.315405

習題EX4-9-2 C++語言

習題EX4-9-2  C++語言

// C++ program to implement Trapezoidal rule
#include<stdio.h>
#include<math.h>

// A sample function whose definite integral's
// approximate value is computed using Trapezoidal
// rule
float y(float x)
{
  // Declaring the function f(x) = x( e^x-1)
    return ( x* (exp(x) -1) );
}

// 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 = 0;
 float xn = 1;

 // Number of grids. Higher value means
 // more accuracy
 int n = 10;

 printf("Value of integral is %6.4f\n",
    trapezoidal(x0, xn, n));
 return 0;
}

設定 

 int n = 10;
結果
Value of integral is 0.5037

設定 

 int n = 100;
結果
Value of integral is 0.5000

習題EX4-9-1 C++語言

習題EX4-9-1 C++語言

// C++ program to implement Trapezoidal rule
#include<stdio.h>
#include<math.h>

// A sample function whose definite integral's
// approximate value is computed using Trapezoidal
// rule
float y(float x)
{
  // Declaring the function f(x) = 0.5x^2
    return ( 0.5*pow(x,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 = 0;
 float xn = 2;

 // Number of grids. Higher value means
 // more accuracy
 int n = 10;

 printf("Value of integral is %6.4f\n",
    trapezoidal(x0, xn, n));
 return 0;
}

輸出結果

Value of integral is 1.3334

範例EX4-8 梯行積分法 雙重積分計算 --- C語言

範例EX4-8 梯行積分法 雙重積分計算 --- C語言

/* ex4-8.c based on Trapezoidal Rule to
 * compute the double integral.
 */
#include <stdio.h>
#include <math.h>
#define F(x,y)  (4*x*y)
#define C(x)    (0)
#define D(x)    (1-pow(x,2))
void gy(int);
double x[50],y[50][50],gg[50],hy[50];

void main()
{
   int i,j,n;
   double a,b,hx,sum=0.0,ts;
   scanf("n=%d a=%lf b=%lf",&n,&a,&b);
   hx=(b-a)/n;
   for(i=0;i<=n;i++)
   {
      x[i]=a+i*hx;
      hy[i]=(D(x[i])-C(x[i]))/n;
      for(j=0;j<=n;j++)
         y[i][j]=C(x[i])+j*hy[i];
   }
   gy(n);
   for(i=1;i<=n-1;i++)
      sum=sum+gg[i];
   ts=(hx/2.0)*(gg[0]+gg[n]+2*sum);
   printf("T%d=%.6lf\n",n,ts);
   return 0;
}

void gy(int n)
{
   int i,j;
   double sum=0.0;
   for(i=0;i<=n;i++)
   {
      for(j=1;j<=n-1;j++)
      {
  sum=sum+F(x[i],y[i][j]);
      }
      gg[i]=(hy[i]/2.0)*(F(x[i],y[i][0])+
     F(x[i],y[i][n])+2*sum);
      sum=0.0;
   }
   return;
}

輸入資料
n=10  a=0.0  b=1.0

輸出結果
T10=0.331650

範例EX4-8 梯行積分法 雙重積分計算-----C++ 與 python 語言

範例EX4-8 梯行積分法 雙重積分計算  C++ 與 python 

  

Trapezoidal Rule for Approximate Value of Definite Integral

In the field of numerical analysis, Trapezoidal rule is used to find the approximation of a definite integral. The basic idea in Trapezoidal rule is to assume the region under the graph of the given function to be a trapezoid and calculate its area.
It follows that:

For more accurate results the domain of the graph is divided into n segments of equal size as shown below:

Grid spacing or segment size h = (b-a) / n.
Therefore, approximate value of the integral can be given by:

// C++ program to implement Trapezoidal rule
#include<stdio.h>
#include<math.h>

// A sample function whose definite integral's
// approximate value is computed using Trapezoidal
// rule
float y(float x)
{
  // Declaring the function f(x) = 2x-4x^3-2x^5
    return ( (2*x)-(4*pow(x,3))+(2*pow(x,5)) );
}

// 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 = 0;
float xn = 1;

// Number of grids. Higher value means
// more accuracy
int n = 10;

printf("Value of integral is %6.4f\n",
trapezoidal(x0, xn, n));
return 0;
}

輸出結果

Value of integral is 0.3316

＝＝＝＝＝＝＝＝＝＝
python 語言

＝＝＝＝＝＝＝＝＝＝

# Python3 code to implement Trapezoidal rule

# A sample function whose definite
# integral's approximate value is
# computed using Trapezoidal rule
def y( x ):

# Declaring the function
# f(x) =  return ( (2*x)-(4*pow(x,3))+(2*pow(x,5)) )
return  ( (2*x)-(4*pow(x,3))+(2*pow(x,5)) )

# Function to evalute the value of integral
def trapezoidal (a, b, n):

# Grid spacing
h = (b - a) / n

# Computing sum of first and last terms
# in above formula
s = (y(a) + y(b))

# Adding middle terms in above formula
i = 1
while i < n:
s += 2 * y(a + i * h)
i += 1

# h/2 indicates (b-a)/2n.
# Multiplying h/2 with s.
return ((h / 2) * s)

# Driver code to test above function
# Range of definite integral
x0 = 0
xn = 1

# Number of grids. Higher value means
# more accuracy
n = 10
print ("Value of integral is ",
"%.4f"%trapezoidal(x0, xn, n))



# This code is contributed by "Sharad_Bhardwaj". 

輸出結果

Value of integral is 0.3316