C 中的梯形黎曼和
Trapezoidal Riemann Sum in C
我一直在尝试用黎曼求和来近似计算 C 中的积分。在我下面的代码中,我尝试用梯形方式和矩形方式来近似(显然,梯形方式应该更好) .
我试着在纸上为此做了一个算法,我得到了以下结果:
注:N 是矩形(或梯形)的数量,dx 是使用 a、b 和 N(dx = (b-a)/N)计算的。 f(x) = x^2
矩形法:
<a href="http://www.codecogs.com/eqnedit.php?latex=\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;f(a&space;+&space;(i-1)dx)dx" target="_blank"><img src="http://latex.codecogs.com/png.latex?\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;f(a&space;+&space;(i-1)dx)dx" title="\int_a^b x^2 dx \approx \sum_{i=1}^N f(a + (i-1)dx)dx" /></a>
梯形法:
<a href="http://www.codecogs.com/eqnedit.php?latex=\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;[f(a&space;+&space;(i-1)dx)&space;+&space;f(a&space;+&space;i\cdot&space;dx)]dx" target="_blank"><img src="http://latex.codecogs.com/png.latex?\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;[f(a&space;+&space;(i-1)dx)&space;+&space;f(a&space;+&space;i\cdot&space;dx)]dx" title="\int_a^b x^2 dx \approx \sum_{i=1}^N [f(a + (i-1)dx) + f(a + i\cdot dx)]dx" /></a>
代码(在下面的代码中,f(x)=x^2 和 F(x) 是反导数 (x^3/3):
int main() {
int no_of_rects;
double a, b;
printf("Number of subdivisions = ");
scanf("%d", &no_of_rects);
printf("a = ");
scanf("%lf", &a);
printf("b = ");
scanf("%lf", &b);
double dx = (b-a)/no_of_rects;
double rectangular_riemann_sum = 0;
int i;
for (i=1;i<=no_of_rects;i++) {
rectangular_riemann_sum += (f(a + (i-1)*dx)*dx);
}
double trapezoidal_riemann_sum = 0;
int j;
for (j=1;j<=no_of_rects;j++) {
trapezoidal_riemann_sum += (1/2)*(dx)*(f(a + (j-1)*dx) + f(a + j*dx));
printf("trapezoidal_riemann_sum: %lf\n", trapezoidal_riemann_sum);
}
double exact_integral = F(b) - F(a);
double rect_error = exact_integral - rectangular_riemann_sum;
double trap_error = exact_integral - trapezoidal_riemann_sum;
printf("\n\nExact Integral: %lf", exact_integral);
printf("\nRectangular Riemann Sum: %lf", rectangular_riemann_sum);
printf("\nTrapezoidal Riemann Sum: %lf", trapezoidal_riemann_sum);
printf("\n\nRectangular Error: %lf", rect_error);
printf("\nTrapezoidal Error: %lf\n", trap_error);
return 0;
}
其中:
double f(double x) {
return x*x;
}
double F(double x) {
return x*x*x/3;
}
我已经包含了 math 和 stdio 头文件。发生的事情是矩形黎曼和没问题,但梯形黎曼和不知为何总是0。
有什么问题?它是我的公式中的东西吗?还是我的代码?
(顺便说一句,我是C的新手)
提前致谢。
在此声明中:
trapezoidal_riemann_sum += (1/2)*(dx)*(f(a + (j-1)*dx) + f(a + j*dx));
1/2
== 零,所以整个语句为零。至少将分子或分母更改为双精度形式以获得双精度值。即 1/2.0
或 1.0/2
或 1.0/2.0
都可以。
我一直在尝试用黎曼求和来近似计算 C 中的积分。在我下面的代码中,我尝试用梯形方式和矩形方式来近似(显然,梯形方式应该更好) .
我试着在纸上为此做了一个算法,我得到了以下结果: 注:N 是矩形(或梯形)的数量,dx 是使用 a、b 和 N(dx = (b-a)/N)计算的。 f(x) = x^2
矩形法:
<a href="http://www.codecogs.com/eqnedit.php?latex=\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;f(a&space;+&space;(i-1)dx)dx" target="_blank"><img src="http://latex.codecogs.com/png.latex?\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;f(a&space;+&space;(i-1)dx)dx" title="\int_a^b x^2 dx \approx \sum_{i=1}^N f(a + (i-1)dx)dx" /></a>
梯形法:
<a href="http://www.codecogs.com/eqnedit.php?latex=\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;[f(a&space;+&space;(i-1)dx)&space;+&space;f(a&space;+&space;i\cdot&space;dx)]dx" target="_blank"><img src="http://latex.codecogs.com/png.latex?\int_a^b&space;x^2&space;dx&space;\approx&space;\sum_{i=1}^N&space;[f(a&space;+&space;(i-1)dx)&space;+&space;f(a&space;+&space;i\cdot&space;dx)]dx" title="\int_a^b x^2 dx \approx \sum_{i=1}^N [f(a + (i-1)dx) + f(a + i\cdot dx)]dx" /></a>
代码(在下面的代码中,f(x)=x^2 和 F(x) 是反导数 (x^3/3):
int main() {
int no_of_rects;
double a, b;
printf("Number of subdivisions = ");
scanf("%d", &no_of_rects);
printf("a = ");
scanf("%lf", &a);
printf("b = ");
scanf("%lf", &b);
double dx = (b-a)/no_of_rects;
double rectangular_riemann_sum = 0;
int i;
for (i=1;i<=no_of_rects;i++) {
rectangular_riemann_sum += (f(a + (i-1)*dx)*dx);
}
double trapezoidal_riemann_sum = 0;
int j;
for (j=1;j<=no_of_rects;j++) {
trapezoidal_riemann_sum += (1/2)*(dx)*(f(a + (j-1)*dx) + f(a + j*dx));
printf("trapezoidal_riemann_sum: %lf\n", trapezoidal_riemann_sum);
}
double exact_integral = F(b) - F(a);
double rect_error = exact_integral - rectangular_riemann_sum;
double trap_error = exact_integral - trapezoidal_riemann_sum;
printf("\n\nExact Integral: %lf", exact_integral);
printf("\nRectangular Riemann Sum: %lf", rectangular_riemann_sum);
printf("\nTrapezoidal Riemann Sum: %lf", trapezoidal_riemann_sum);
printf("\n\nRectangular Error: %lf", rect_error);
printf("\nTrapezoidal Error: %lf\n", trap_error);
return 0;
}
其中:
double f(double x) {
return x*x;
}
double F(double x) {
return x*x*x/3;
}
我已经包含了 math 和 stdio 头文件。发生的事情是矩形黎曼和没问题,但梯形黎曼和不知为何总是0。
有什么问题?它是我的公式中的东西吗?还是我的代码? (顺便说一句,我是C的新手)
提前致谢。
在此声明中:
trapezoidal_riemann_sum += (1/2)*(dx)*(f(a + (j-1)*dx) + f(a + j*dx));
1/2
== 零,所以整个语句为零。至少将分子或分母更改为双精度形式以获得双精度值。即 1/2.0
或 1.0/2
或 1.0/2.0
都可以。