Вычисление интеграла составной квадратурой Гаусса с 2 узлами

[aristarty]

Добрый день!

Разбираюсь с заданием - вычислить с помощью квадратурной формулы Гаусса с двумя узлами интеграл Френеля на отрезке от 0 до 1.5 при разных разбиениях этого интервала.

$C(x) = \int_0^{x}{\cos (\frac{\pi t^2}{2})dt}$ - этот интеграл при х=1.5 я также считаю с помощью разложения в ряд Тейлора

Я задаю программе количество разбиений N, строю разбиение отрезка [0,1.5] на N частей - это Z-массив.

Далее вычисляю с помощью составной квадратурной формулы интеграл следующим образом:

$\int_c^{d}{\varphi(t)dt}= \sum_{i=1}^{N}{\int_{z_{i-1}}^{z_i}{\varphi(t)dt}}\approx \sum_{i=1}^{N}{S_i(\varphi)}$

где $S_i(\varphi) = \frac{h_N}{2}[\varphi(z_{i-1}+\frac{h_N}{2}(1-\frac{1}{\sqrt{3}}))+\varphi(z_{i-1}+\frac{h_N}{2}(1+\frac{1}{\sqrt{3}}))]$

В ходе расчетов я анализирую модуль расхождения значений интеграла, посчитанного с помощью ряда Тейлора и с помощью квадратурной формулы.

$\left | C(x)- \sum_{i=1}^{N}{S_i(\varphi)} \right |$

Преподаватель как-то сказал, для квадратуры Гаусса c 2 узлами для разных разбиений N справедливо $\frac{\left | C(x)- \sum_{i=1}^{N}{S_i(\varphi)} \right |}{{h_N}^{4}} = \operatorname{const}$

Это своеобразная проверка программы на вшивость. К сожалению, не могу найти этому подтверждение. И долго пытаюсь понять, где в ходе рассчетов я поступаю неправильно? Поскольку $\frac{\left | C(x)- \sum_{i=1}^{N}{S_i(\varphi)} \right |}{{h_N}^{4}} $ принимает разные значения

Например

-------------
x=1,5
N=16
Гаусс = 0,5368527
Интеграл = 0,4452612
Разница = 0,09159157
Разница/h^4 = 1185,688
------------
x=1,5
N=32
Гаусс = 0,4901777
Интеграл = 0,4452612
Разница = 0,04491657
Разница/h = 9303,396
-------------
x=1,5
N=100
Гаусс = 0,4452625
Интеграл = 0,4452612
Разница = 1,311302E-06
Разница/h = 25,90227
-------------

Код программы прилагаю:

//Rextester.Program.Main is the entry point for your code. Don't change it.
//Microsoft (R) Visual C# Compiler version 2.9.0.63208 (958f2354)

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text.RegularExpressions;

namespace Rextester
{
    public class Program
    {
        
        public static void getGrid(float a, float b, float step, ref List<float> X)
        {
            float x = a;
            

            while (x < b)
            {
                X.Add(x);
                x += step;

            }

        }
        
        public static float Cos_count(float val)
        {
            return (float)Math.Cos(Math.PI * val * val / 2);   //вычисление подынтегральной функции интеграла Френеля
        }
        
        public static float S(int N, float x) {   //расчет квадратуры Гаусса
        
            List<float> Z = new List<float>();
            float h = x / N;
            Program.getGrid(0, x, h, ref Z);  //получение разбиения

            float valS = 0;
            float S_i = 0;
            for (int i=1; i<Z.Count(); i++)
            {
                S_i = (h / 2) * (Program.Cos_count((float)(Z[i-1]+(h/2)*(1-1/Math.Sqrt(3))))+ Program.Cos_count((float)(Z[i - 1] + (h / 2) * (1 + 1 / Math.Sqrt(3)))));
                valS += S_i;
            }
            return valS;
        }
        
        public static float Integral_counter(float b, float ep)  //расчет интеграла по Тейлору
        { 
        }
        
        
        public static float GaussCounter (int n, float x){
            float eps = 1e-4F;
       
            float S_n = S(n,x);
            //float S_2n = S(2 * n,x);

              //  while (Math.Abs(S_2n - S_n) > eps)
                //{
                  //  S_n = S_2n;
                   // n *= 2;
                   // S_2n = S( 2 * n,x);
               // }
            return S_n;
        }
        

        public static void Main(string[] args)
        {
            int n=100;
            float x=1.5F;
   
            float step=(1.5F)/n;
                      
            Console.WriteLine("-------------");
            Console.WriteLine("x="+x);          
            Console.WriteLine("N="+n);          
                
            Console.WriteLine("Гаусс = "+GaussCounter(n,x));
            Console.WriteLine("Интеграл = "+Integral_counter(x,1e-6F));
                Console.WriteLine("Разница = "+Math.Abs(Integral_counter(x,1e-6F)-GaussCounter(n,x)));
                Console.WriteLine("Разница/h = "+Math.Abs(Integral_counter(x,1e-6F)-GaussCounter(n,x))/((x/n)*(x/n)*(x/n)*(x/n)));
            Console.WriteLine("-------------");
           
            
            
        }
    }
}

И вообще корректна ли "проверка" программы с помощью отслеживания константности расхождения, деленного на шаг в четвертой степени в данном случае?

Спасибо заранее всем за любые адекватные идеи.

[TOTAL]

Цитата:

for (int i=1; i<Z.Count(); i++)
{
S_i = (h / 2) * (Program.Cos_count((float)(Z[i-1]+(h/2)*(1-1/Math.Sqrt(3))))+ Program.Cos_count((float)(Z[i - 1] + (h / 2) * (1 + 1 / Math.Sqrt(3)))));

Уберите подпрограмму с сеткой, массив для сетки, величины сделайте 8 байтовыми и замените кусок на

  for (int i=0; i<N; i++)
            {
                      z0=h*(i+1./2);   dd=h/(2*Math.Sqrt(3.))
              S_i = h / 2 * (Program.Cos_count(z0-dd)+ Program.Cos_count(z0+dd));

[aristarty]

TOTAL
Спасибо большое за предложенные исправления, расхождение стало гораздо меньше (~ до 9 знака совпадают). Но, к сожалению, константности расхождения, деленного на четвертую степень шага разбиения, по-прежнему нет, не подскажете в какую сторону копать по этому вопросу?

Прим. для тех же N, что и в первом сообщении, получается следующее:

x=1,5
N=16
Гаусс = 0,445261146183047
Интеграл = 0,445261174682689
Разница = 2,84996415911287E-08
Разница/h^4 = 0,0003689387676674
---------

x=1,5
N=32
Гаусс = 0,44526117474882
Интеграл = 0,445261174682689
Разница = 6,6131045084461E-11
Разница/h^4 = 1,36974670084906E-05
-----------

x=1,5
N=100
Гаусс = 0,445261176028063
Интеграл = 0,445261174682689
Разница = 1,34537442297855E-09
Разница/h^4 = 0,0265752972440208
-------------

[TOTAL]

aristarty в сообщении #1515626 писал(а):

TOTAL
Спасибо большое за предложенные исправления, расхождение стало гораздо меньше (~ до 9 знака совпадают). Но, к сожалению, константности расхождения, деленного на четвертую степень шага разбиения, по-прежнему нет, не подскажете в какую сторону копать по этому вопросу?

Стремление к константе при увеличении числа разбиений наблюдалось бы, если бы вычисления были точными, а не машинными. Поэтому делайте расчеты при маленьких $N=5,10,15,20,25,30$. То есть когда погрешность квадратурной формулы больше погрешности машинных вычислений.

[ewert]

TOTAL в сообщении #1515673 писал(а):

То есть когда погрешность квадратурной формулы больше погрешности машинных вычислений.

При $n=100$ за счёт округлений мы теряем, грубо говоря, один десятичный порядок, что на фоне наблюдаемой нерегулярности совершенно несущественно.

[aristarty]

Я считаю значение функции по ряду Тейлора с точностью 1е-6, уменьшение точности до 1е-16 позволило избежать проблем. Спасибо всем за ответы!