0

我正在使用二分法在从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;
}
4

0 回答 0