Интегралы

// z2-5.cpp : Defines the entry point for the console application.
//
 
#include "stdafx.h"
 
 
int _tmain(int argc, _TCHAR* argv[])
{
    return 0;
}
 
/*Вычисление определенного интеграла функции sqrt (1+cos^2(x))*/
 
/*ЭТАЛОН 3.8202*/
 
#include <stdio.h>
#include <math.h>
#define M_PI 3.14159265358979323846 //В Visual'e M_PI не пашет, к сожалению, приходится выкручиваться
 
double f_value (double x); //Значение подинтегрального выражения
 
double pryam (double lef, double rig, int N); //Вычисление методом прямоугольника
double trap (double lef, double rig, int N); //Вычисление методом трапеций
double simp (double lef, double rig, int N);//Вычисление методом Симпсона
 
double PryamAcc (double lef, double rig, double acc, int *pNp); //Вычисление методом прямоугольника с заданной точностью
double TrapAcc (double lef, double rig, double acc, int *pNt); //Вычисление методом трапеций с заданной точностью
double SimpAcc (double lef, double rig, double acc, int *pNs); //Вычисление методом Симпсона с заданной точностью
 
int main (void)
{
    double lef = 0.0, rig = M_PI, eps; //Где lef(left) - нижняя граница интеграла, rig(right) - правая, eps(epsilon) - точность
    double val_p, val_t,val_s;  //Значения каждого метода с заданной точностью
    int Np, Nt, Ns, i = 1;      //Числа итераций и порядок точности
    
    printf ("E(acc)|Rectangles| N | Trapezoid| N | Simpson  | N\n"); 
 
    for (eps = 0.1; eps >= 0.00001; eps /= 10, i++)
    {
        Np = 1; Nt = 1; Ns = 1;
 
        val_p = PryamAcc (lef, rig, eps, &Np);
        val_t = TrapAcc (lef, rig, eps, &Nt);
        val_s = SimpAcc (lef, rig, eps, &Ns);
 
        printf("10e-%d | %lf | %d | %lf | %d | %lf | %d \n", i, val_p, Np, val_t, Nt, val_s, Ns);
    }
 
    scanf ("%lf", &lef); //Отладка
    
    return 0;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
double PryamAcc (double lef, double rig, double acc, int *pNp) //С заданной точностью
{
    double res, res1, di; //Значение интеграла, уточненное значение интеграла, разность интегралов в каждом цикле 
    
    res = pryam(lef, rig, (*pNp));
    do
    {
        (*pNp) += 1;    //Для красоты
        //(*pNp) *= 2;  //Обычно, т.к точность пропорциональна шагу
        res1 = pryam(lef, rig, (*pNp));
        di = fabs(res - res1);
        res = res1;
    }
    while(di > acc);
    return res1;
}
 
double TrapAcc (double lef, double rig, double acc, int *pNt) //С заданной точностью
{
    double res, res1, di; //Значение интеграла, уточненное значение интеграла, разность интегралов в каждом цикле 
    
    res = trap(lef, rig, (*pNt));
    do
    {
        (*pNt) += 1;    //Для красоты
        //(*pNp) *= 2;  //Обычно, т.к точность пропорциональна шагу
        res1 = trap(lef, rig, (*pNt));
        di = fabs(res - res1);
        res = res1;
    }
    while(di > acc);
    return res1;
}
 
double SimpAcc (double lef, double rig, double acc, int *pNs) //С заданной точностью
{
    double res, res1, di; //Значение интеграла, уточненное значение интеграла, разность интегралов в каждом цикле 
    
    res = simp(lef, rig, (*pNs));
    do
    {
        //(*pNs) += 1;  //Для красоты
        (*pNs) *= 2;    //Обычно, т.к точность пропорциональна шагу
        res1 = simp(lef, rig, (*pNs));
        di = fabs(res - res1);
        res = res1;
    }
    while(di > acc);
    return res1;
}
 
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
double f_value (double x) //Вычисление функции
{
        return sqrt(1.0+cos(x)*cos(x)); //Значение функции
        //return x;
        //return cos(x);
}
 
double pryam (double lef, double rig, int N) //Метод прямоугольников
{
    double dx, x, summ = 0.0; //Шаг сетки; аргумент; вычисленное (возвращаемое) значение интеграла
    dx = (rig-lef)/N;
    for (x = lef; x < rig - dx*0.9; x += dx)
    {
        summ += f_value(x);
    }
    return summ*dx;
}
 
double trap (double lef, double rig, int N) //Метод  трапеции
{
    double dx, x, summ = 0.0; //Шаг сетки; аргумент; вычисленное (возвращаемое) значение интеграла
    dx = (rig-lef)/N;
    summ = 0.5*f_value(lef)+0.5*f_value(rig); //(f(x0)-f(xN))/2
    for (x = lef + dx; x < rig - dx*0.9; x += dx)
    {
        summ += f_value(x);
    }
    return summ*dx;
}
 
double simp (double lef, double rig, int N) //Метод Симпсона
{
    double dx, x, summ1 = 0.0, summ2 = 0.0; //Шаг сетки; аргумент; вычисленное (возвращаемое) значение интеграла
    dx = (rig-lef)/N;
    for (x = lef + 2*dx; x < rig-dx*0.9; x += 2*dx) //x с четным индексом (при lef = x0)
    {
        summ1 += f_value(x);
    }
    for (x = lef + dx; x < rig - dx*0.9; x += 2*dx) //x с нечетным индексом (при lef = x0)
    {
        summ2 += f_value(x);
    }
    return (f_value(lef) + f_value(rig) + 2*summ1 + 4*summ2)*dx/3;
}