Вычислить интеграл по формуле Симпсона - C (СИ)

Формулировка задачи:

Требуется вычислить интеграл по формуле Симпсона. (с экономией вычисления) Наработки:

#include <stdio.h>
#include <math.h>
double Simpson (double a, double b, double eps, double (*w)(double));
double f(double);
main ()
{
    double s,a,b,eps;
    printf ("\n Zadaite interval i tochnost':");
    scanf("%lf%lf%lf",&a,&b,&eps);
    s=Simpson(a,b,eps,f);
    printf("\n Integral ot a=%3.2if do b=%3.2if raven %lf",a,b,s);
}
double f(double x)
{
    return sqrt(1+pow(sin(x),3));
}
double Simpson (double a, double b, double eps, double(*w)(double))
{
    double h,sn,s1,s2,s4,s;
    int n,i;
    n=2;
    s1=w(a)+w(b);
    s2=0;
    s4=w(.5*(a+b));
    s=1e25;
    do {
        sn=2; n*=2; s2+=s4;
        h=(b-a)/n;
        s4=0;
        for (i=0;i<n/2;i++)
            s4+=w(a+(2*i+1)*h);
        s=h/3*(s1+2*s2+4*s4);
    } while (fabs(s-sn)>eps);
    return sn;
    }

Решение задачи: «Вычислить интеграл по формуле Симпсона»

Листинг программы

double Simpson(double a, double b, double dx, double eps)
{
    double x = a, I1, I2 = 0, s;
 
    do{
        I1 = I2;
        s = (f(a) - f(b))/2;
        x = a + dx;
 
        while(x < b){
            s += 2 * f(x) + f(x + dx);
            x += 2 * dx;}
        I2 = 2 * dx * s / 3;
        dx /= 2.0;
    }while(fabs(I1 - I2) > eps);
 
    return I2;
}

Объяснение кода листинга программы

Входные параметры функции: a, b, dx, eps.
Переменная x инициализируется значением a.
Переменная I1 инициализируется значением 0.
Переменная I2 инициализируется значением 0.
Переменная s инициализируется значением (f(a) - f(b))/2.
Первый элемент списка (a, b, dx, eps) используется для вычисления значения переменной s.
Используя цикл do-while, происходит последовательное улучшение приближения до тех пор, пока изменение интеграла не станет меньше заданной точности eps.
Цикл do-while состоит из двух циклов: внешний цикл do-while и внутренний цикл while.
Внутренний цикл while используется для вычисления среднего значения функции на каждом элементе списка (x, dx).
Значение переменной s обновляется на каждом шаге внутреннего цикла while путем добавления суммы двух значений функции и значения функции на следующем шаге.
Переменная I2 обновляется на каждом шаге внешнего цикла do-while путем умножения значения переменной s на 2 * dx / 3.
Значение переменной dx уменьшается в два раза на каждом шаге внешнего цикла do-while.
Переменная x обновляется на каждом шаге внутреннего цикла while путем добавления значения dx к текущему значению x.
Значение переменной dx обновляется на каждом шаге внешнего цикла do-while путем деления на 2.0.
Значение переменной I1 обновляется на каждом шаге внешнего цикла do-while путем присвоения значения переменной I2.
Выходные параметры функции: возвращаемое значение I2.
Функция Simpson завершается, когда разница между последними двумя приближениями I1 и I2 становится меньше заданной точности eps.
Значение переменной I2 возвращается как результат функции.
Значение переменной I2 является приближенным значением интеграла по формуле Симпсона.
Значение переменной eps используется для контроля точности приближения.