我正在使用二分法在从70*10^9到250*10^9的域中找到函数的根,但输出始终是上限,即250*10^9。函数是定积分,不知道哪里做错了。提前致谢。
#include <stdio.h>
#include <math.h>
double I = 0.0225, W = 50000, L = 25, Y = 12000;
double x, E;
typedef double(*DfD)(double);
double g(double), h(double);
double pow(double a, double b);
double g(double x){
return W*x*(x-L/2)-(Y*pow(x,3))/2;
}
double h(double x){
return Y*pow(x,3)/2;
}
double midpoint_int(DfD f, double x0, double x1, int n) {
int i;
double x, dx, sum = 0.0;
dx = (x1-x0)/n;
for (i = 0, x = x0 + dx/2; i < n; i ++, x +=dx)
sum += f(x);
return sum*dx;
}
double deflection (DfD f, double int1, double int2){
int1 = midpoint_int(g, L/2, L, 1000);
int2 = midpoint_int(h, 0, L/2, 1000);
return (int1 - int2)/(E*I)+0.5;
}
double bisection(DfD f, double a0, double a1, double tol ){
double middle;
for (;;){
middle = (a0 + a1)/2.0;
if (fabs(middle - a0) < tol)
return middle;
else if (f(middle) * f(a0) < 0.0)
a1 = middle;
else
a0 = middle;
}
}
int main (void){
double E;
E = bisection (deflection, 70*pow(10,9), 250*pow(10,9), 0.001);
printf("The optimal elastic modulus is %fPa.\n", E);
return;
}