метод симпсона и интегралы в си

[Михаил 125]

Снимок экрана 2022-12-11 в 19.46.46.png
double Simpson(double, double, double, double(*) (double));
double f(double);
int main()
{
double s,a,b,eps;
printf ("\n Задайте пределы интегрирования и точность:");
scanf ("%lf%lf%lf",&a,&b,&eps);
s=Simpson (a,b,eps,f);
printf ("\n Интеграл от a=%3.21f до b=%3.21f равен %lf",a,b,s);
}
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=s;n*=2; s2+=s4;
h=(b-a)/n;
s4=0;
for (i=0;i<n;i++)
s4+=w(a+(2*i+1)*h);
s=h/3*(s1+2*s2+4*s4);
}
while (fabs(s-sn)>eps);
return sn;
}
double f (double x)
{
return 1/sqrt((1+pow(sin(x),3)));
}
При вводе значений 1 2 и точности 0.1 выдает значение 1.4.... если же решать пример через онлайн калькулятор, то выдает 0.7...
Проверьте, пожалуйста, может я не те значения ввожу... и подскажите как можно посчитать интеграл с разбиением на 4 и на 8 частей

[type_Oleg]

Вот, взял отсюда https://ru.stackoverflow.com/questio...BE%D0%BD%D0%B0 , немного изменил.

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

double fun(double x)
{
	return 1/sqrt(1+pow(sin(x),3));
}
double Simpson(double a, double b, double eps, double (*f)(double))
{
    double N = 4.;
    double sum1 = 0., sum2 = 0.;
    double x, h;
    double index1, index2;
    unsigned int i;
    double func;
    for(;;) {  //  можно задать макс число циклов на всякий случай
        h = (b-a)/N;
        sum1 = sum2 = -(f(a)+f(b));
        index1 = index2 = 4.;
        i = 0;
        for(x = a; x <= b+h/2.; x+=h) {
            sum1 +=  (index1 = 6.-index1)*(func=f(x));
            if (i++%2==0) {
                sum2 += (index2 = 6.-index2)*func;
            };
        };
        if (fabs((sum1-2.*sum2)/sum1)<=eps) return (h*sum1)/3.+h*(2.*sum2-sum1)/45.; else N*=2.;
    };
    return (h*sum1)/3.;
};

int main(int argc, char **argv)
{
	double s;
	s = Simpson(1.0,2.0,0.00001,fun);
	printf (" %lf",s);
	return 0;
}

Выдает 0.73033. Проверял, это правильный результат.

[Михаил 125]

type_Oleg, спасибо, а как можно посчитать интеграл с разбиениями на 4 и на 8 частей ?

[type_Oleg]

Сообщение от Михаил 125

type_Oleg, спасибо, а как можно посчитать интеграл с разбиениями на 4 и на 8 частей ?

Я сам не знаю, какой метод я вам предложил. Просто скопировал, поправил, проверил.
Может быть и есть на 4 части -

  double N = 4.;
....
  h = (b-a)/N;

Уж сами разбирайтесь.

PS
Кстати, вроде зря

  h = (b-a)/N;

внутри цикла. Вроде a,b,N не меняются в цикле, тогда зачем одно и то же h вычислять в цикле.