#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]);

}

}