《計算方法》 李曉紅等 第六章 數值積分 例題 程式碼
#include "stdafx.h"
#include <math.h>
// 數值積分
// Newton-Cotes積分公式
void
NewtonCotesIntegral()
{
printf
(
"牛頓-科茨積分公式:\n"
);
int
n;
double
h;
int
i;
double
a=0,b=0.5;
double
I;
double
C[5],f[5],x[5];
for
(n=1;n<=4;n*=2)
{
h = (b-a)/n;
for
(i=0;i<=n;i++) {
x[i] = a+i*h;
// 求積節點
f[i] = 1./(1+x[i]*x[i]);
// 節點函式值
}
// Cotes係數
if
(n==1)
C[0] = C[1] = 0.5;
else
if
(n==2)
{
C[0] = C[2] = 1./6;
C[1] = 4./6;
}
else
if
(n==4)
{
C[0] = C[4] = 7./90;
C[1] = C[3] = 16./45;
C[2] = 2./15;
}
I = 0.;
for
(i=0;i<=n;i++)
I += C[i]*f[i];
I *= (b-a);
printf
(
"n=%d,I=%.9f \n"
,n,I);
}
printf
(
"真值I=%.9f\n"
,
atan
(0.5));
}
void
FuHuaTiXingIntegral()
{
printf
(
"復化梯形求積:\n"
);
int
n[4] = {1,2,4,5};
double
a=0,b=0.5;
int
i,j;
int
nn;
double
h;
double
x[6],f[6];
double
I;
for
(i=0;i<4;i++) {
nn = n[i];
h = (b-a)/nn;
for
(j=0;j<=nn;j++) {
x[j] = a+j*h;
f[j] = 1./(1+x[j]*x[j]);
}
// PS(x,nn+1,"x");
// PS(f,nn+1,"f");
I = 0;
for
(j=1;j<=nn;j++) {
I += (f[j-1]+f[j])*h/2;
}
printf
(
"n=%d,I=%.9f\n"
,nn,I);
}
}
void
FuHuaSimpsonIntegral()
{
printf
(
"復化辛普生求積:\n"
);
int
m,n;
double
h;
double
a=0,b=0.5;
double
x[7],f[7];
int
i;
double
I;
for
(m=1;m<=3;m++) {
n=2*m;
h = (b-a)/n;
for
(i=0;i<=n;i++) {
x[i] = a+i*h;
f[i] = 1./(1+x[i]*x[i]);
}
I = 0;
for
(i=0;i<m;i++) {
I += (f[2*i]+4*f[2*i+1]+f[2*i+2]);
}
I = I*h/3;
printf
(
"m=%d,I=%.9f\n"
,m,I);
}
}
// 復化科茨、辛普森、梯形積分
void
FuHuaCotesIntegral()
{
printf
(
"復化科茨、辛普森、梯形積分:\n"
);
int
m = 4;
int
n = 2*m;
int
i,j;
double
T8,S8,C8;
double
x[9],f[9];
double
h;
double
a=0,b=1;
double
C[5];
h = (b-a)/n;
for
(j=0;j<=n;j++) {
x[j] = a+j*h;
if
(j==0) f[j]=1;
else
f[j] =
sin
(x[j])/x[j];
}
// 梯形T8
T8 = 0.;
for
(j=1;j<=n;j++) {
T8 += (f[j-1]+f[j]);
}
T8 *= h/2;
printf
(
"n=%d,T8=%.9f\n"
,n,T8);
// 辛普森
S8 = 0.;
for
(i=0;i<m;i++) {
S8 += (f[2*i]+4*f[2*i+1]+f[2*i+2]);
}
S8 = S8*h/3;
printf
(
"m=%d,S8=%.9f\n"
,m,S8);
// 科茨(4)
C8 = 0.;
C[0] = C[4] = 7./90;
C[1] = C[3] = 16./45;
C[2] = 2./15;
for
(i=0;i<n/4;i++)
{
for
(j=0;j<=4;j++)
{
C8 += (C[j]*f[4*i+j]);
}
}
C8 *= (b-a)/(n/4);
printf
(
"m=%d,C8=%.9f\n"
,8,C8);
}
void
ChangeStepIntegral()
{
printf
(
"變步長數值積分: 梯形法\n"
);
int
n = 1;
int
k = 0;
double
T0,T1,RT8;
double
a= 0,b=1.;
double
h;
int
i;
double
x0,x1,f0,f1;
for
(;;)
{
h=(b-a)/n;
x0 =a; f0 = 1.;
T1 = 0.;
for
(i=1;i<=n;i++)
{
x1=a+i*h;
f1=
sin
(x1)/x1;
T1 += (f0+f1);
f0 = f1;
}
T1 *= (h/2.);
printf
(
"T%d=%.9f\n"
,n,T1);
if
(n==8) RT8=(T1-T0)/3;
k++;
n*=2;
if
(k>10)
break
;
else
T0 = T1;
}
printf
(
"R[f,T8]=%.9f\n"
,RT8);
}
void
RombergIntegral()
{
printf
(
"龍貝格積分: \n"
);
// Ti: 1,2,4,8,16,32,64,...
// Si: 1,2,4, 8,16,32,...
// Ci: 1,2, 4,8,16,...
// Ri: 1, 2,4,8,...
double
f;
double
T[16+1],S[8+1],C[4+1],R[2+1];
int
k,i;
double
h;
double
a=0,b=1.;
int
n=16;
for
(k=0;k<=16;k++) T[k] = 0;
for
(k=0;k<=8;k++) S[k] = 0;
for
(k=0;k<=4;k++) C[k] = 0;
for
(k=0;k<=2;k++) R[k] = 0;
k = 1;
for
(;;)
{
T[k] = 0;
h = (b-a)/k;
for
(i=0;i<=k;i++) {
f = 4./(1.+(a+i*h)*(a+i*h));
if
(i==0||i==k)
T[k] += f;
else
T[k] += 2*f;
}
T[k] *= (h/2);
if
(k%2==0)
S[k/2] = (4*T[k]-T[k/2])/(4-1);
if
(k%4==0)
C[k/4] = (16*S[k/2]-S[k/4])/(16-1);
if
(k%8==0)
{
R[k/8] = (64*C[k/4]-C[k/8])/(64-1);
if
(k/8>1&&
abs
(R[k/8]-R[k/8/2])<1e-5)
break
;
}
k*=2;
}
PS(&T[1],n,
"T"
);
PS(&S[1],n/2,
"S"
);
PS(&C[1],n/4,
"C"
);
PS(&R[1],n/8,
"R"
);
printf
(
"誤差 = %.2e\n"
,
abs
(R[2]-R[1]));
}
// 高斯求積公式
void
GaussLegendreIntegral()
{
printf
(
"高斯-勒讓德求積公式:\n"
);
double
A1[2],A2[3];
double
x1[2],x2[3];
double
a=0,b=1;
int
n,i;
double
G[2];
double
*x,*A;
A1[0]=A1[1]=1;
x1[0]=-1./
sqrt
(3.); x1[1]=-x1[0];
x2[0]= -
sqrt
(0.6);x2[1]=0;x2[2]=-x2[0];
A2[0]=A2[2]=0.5555556;
A2[1]=0.8888889;
for
(n=1;n<=2;n++) {
G[n-1] = 0;
if
(n==1)x=x1;
else
x=x2;
if
(n==1)A=A1;
else
A=A2;
for
(i=0;i<=n;i++) {
G[n-1]+=A[i]*
exp
(x[i]*(b-a)/2+(b+a)/2);
}
G[n-1] *= (b-a)/2;
printf
(
"g%d=%.12f\n"
,n+1,G[n-1]);
}
}
void
FuHuaGaussLegendreIntegral()
{
printf
(
"復化高斯-勒讓德積分:\n"
);
double
A[2];
double
x[2];
double
a=0,b=1;
int
n,i,m;
double
G[3];
double
h;
double
x0,x1;
A[0]=A[1]=1;
x[0]=-1./
sqrt
(3.); x[1]=-x[0];
for
(m=1;m<=3;m++)
// Gm, 區間等分個數
{
h = (b-a)/m;
G[m-1] = 0;
for
(n=1;n<=m;n++)
// 區間
{
x0 = a+(n-1)*h;
x1 = a+n*h;
for
(i=0;i<=1;i++)
// 節點
{
G[m-1] += A[i]*
exp
(x[i]*(x1-x0)/2+(x0+x1)/2);
}
}
G[m-1] *= (h/2);
printf
(
"g%d=%.12f\n"
,m,G[m-1]);
}
}