6-8 Shape Roof (50分)

The kind of roof shown in Figure 1 is shaped from plain flat rectangular plastic board in Figure 2.
Figure 1
Figure 2
The transection of the roof is a sine curve which fits the function y(x)=lsin(lx) in centimeters. Given the length of the roof, your task is to calculate the length of the flat board needed to shape the roof.

Format of function:
double Integral(double a, double b, double (*f)(double x, double y, double z), double eps, double l, double t);
where the function Integral returns the value of a definit integral of a given function pointed by f over a given interval from double a to double b. The integrant function f will be defined in the sample program of judge. Other parameters are: double eps, the accuracy of the resulting value of the integral; double l and double t, the altitude and the frequency of the sine curve, respectively.

Note: the length of the flat board which Integral returns must be in meters.

Sample program of judge:

#include <stdio.h>
#include <math.h>

double f0( double x, double l, double t )
{
    return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x));
}

double Integral(double a, double b, double (*f)(double x, double y, double z), 
double eps, double l, double t);

int main()
{
    double a=0.0, b, eps=0.005, l, t;

    scanf("%lf %lf %lf", &l, &b, &t);
    printf("%.2f\n", Integral(a, b, f0, eps, l, t));

    return 0;
}

/* Your function will be put here */

Sample Input:
2 100 1
Sample Output:
1.68

解法

只需要把条件换成达到eps输出就可以了。

需要注意的一点是：必须先算一个周期，再算小于1个周期的部分。

也就是$(0,T),(a+n\cdot T,b)$这两个区间分别算，但是误差要加起来算。$n\cdot (0,T)+(a+n\cdot T,b)$

上代码

double Integral(double a, double b, double (*f)(double x, double y, double z), double eps, double l, double t)
{
    
    double T = 3.1415926535898 / t;
    int n = (b - a) / T;
    double R1[3][10000];
    double R2[3][10000];
    double h1 = T;
    double h2 = b - a - n * T;
    R1[1][1] = h1 / 2 * (f(0, l, t) + f(T, l, t));
    R2[1][1] = h2 / 2 * (f(a + n * T, l, t) + f(b, l, t));
    for (int i = 2; ; i++)
    {
        double sum1 = 0;
        double sum2 = 0;
        for (int k = 1; k <= pow(2, i - 2); k++)
        {
            sum1 += f(0 + (k - 0.5) * h1, l, t);
            sum2 += f(a + n * T + (k - 0.5) * h2, l, t);
        }
        R1[2][1] = 1.0 / 2 * (R1[1][1] + h1 * sum1);
        R2[2][1] = 1.0 / 2 * (R2[1][1] + h2 * sum2);
        for (int j = 2; j <= i; j++)
        {
            R1[2][j] = R1[2][j - 1] + (R1[2][j - 1] - R1[1][j - 1]) / (pow(4, j - 1) - 1);
            R2[2][j] = R2[2][j - 1] + (R2[2][j - 1] - R2[1][j - 1]) / (pow(4, j - 1) - 1);
        }
        if (fabs(R1[1][i - 1] * n + R2[1][i - 1] - R1[2][i] * n - R2[2][i]) < eps)
            return (R1[2][i] * n + R2[2][i]) / 100;
        h1 = h1 / 2;
        h2 = h2 / 2;
        for (int j = 1; j <= i; j++)
        {
            R1[1][j] = R1[2][j];
            R2[1][j] = R2[2][j];
        }
    }
}