cbsw
2013-11-28 15:58:51 +08:00
Simpson和梯形公式计算积分
/* C version */
float f(float x){
return sqrt(x)*log(x);
}
while (N < 100000){ /* Simpson integration*/
h = (b-a)/N;
Sn = h*f(b)/6;
for (int i = 0; i < N; i++){
if (i == N-1)
Sn += h/6 * (4*f(a+(i+1.0/2)*h));
else
Sn += h/6 * (4*f(a+(i+1.0/2)*h) + 2*f(a+(i+1.0)*h));
}
errS = fabs(Sn-rval);
N += 100;
/* fprintf (fp,"%d\t%.9f\n",N,errS); */
}
N=1000;
while (N<100000){ /* Trapezoidal integration */
h = (b-a)/N;
Tn = 0.0;
for (int i = 0; i < N; ++i) {
if (i!=0)
Tn += h/2 * (f(a+i*h) + f(a+(i+1)*h));
}
errT = fabs(Tn-rval)
N += 100;
/* fprintf(fp,"%d\t%.9f\n",N,errT); */
}
这个程序是用来确定舍入误差开始超过积分方法误差时的积分步长，当然这个程序实际中没有多大意义
# python version
for N in n:
h = (b-a)/N
# Sn
Sn = h*f(b)/6
for k in range(0,N,1):
if k==N-1:
Sn += h/6 * (4*f(a+(k+1.0/2)*h))
else:
Sn += h/6 * (4*f(a+(k+1.0/2)*h) + 2*f(a+(k+1)*h))
# Tn
Tn = 0.0
for k in range(0,N,1):
if k!=0:
Tn += h/2 * (f(a+k*h) + f(a+(k+1)*h))
errS.append(abs(Sn-real))
errT.append(abs(Tn-real))
当然这个程序可以优化，但相同的算法，用不同语言、解释器实现，这样比较应该还算比较公平