Динамический метод численного интегрирования Гаусса‑Лежандра на C#: вычисляем узлы и веса «на лету» без переполнения

2026-10-31 в 12:15, admin, рубрики: метод гаусса, метод Гаусса-Лежандра, тестирование, численное интегрирование, численные методы

Теоретическая часть

Я тут не буду давать всю математику по полиномам Лежандра, а дам, только ту часть их теории, которая непосредственно касается использования их при численном интегрировании. Кто хочет узнать о этих полиномах поподробнее, могут найти учебниках и в интернете.

Полиномы Лежандра n степени являются многочленами n‑ой степени следующего вида:

$P_nleft(xright)=frac{1}{2^n}cdotsum_{k=0}^{left[n/2right]}{left(-1right)^kcdot C_k^ncdot C_n^{2cdot n-2cdot k}cdot x^{n-2cdot k}},,,,, (1)$

Где $C_n^k=frac{n!}{k!cdotleft(n-kright)!}$ — число сочетаний из по k (биномиальный коэффициент), а [n] — целая часть числа.

Формула (1) для полиномов Лежандра очень громоздкая и для реализации на компьютере имеет определенные сложности по причине округления и при вычислении коэффициентов полиномов. Если вы будете вычислять коэффициенты для полиномов Лежандра в лоб, то для полиномов выше 15-й степени столкнетесь с численным переполнением. В этом случае надо реализовать альтернативные вычисления. Самый эффективный альтернативный метод является метод вычисления значений полиномов Лежандра по рекурсивной формуле (соотношению Бонне):

$P_{n+1}left(xright)=frac{2cdot n+1}{n+1}cdot xcdot P_nleft(xright)-frac{n}{n+1}cdot P_{n-1}left(xright)$

Производная полинома Лежандра n‑й степени вычисляется по следующей формуле:

$P_n^primeleft(xright)=frac{n}{1-x^2}cdotleft(P_{n-1}-xcdot P_nleft(xright)right)$

Так как с помощью данной формулы нельзя определить значения производных при x=-1 и x=1 (знаменатель обращается в ноль), то производные полиномов Лежандра вычисляются по следующим формулам:

При x=-1:

$P_n^primeleft(-1right)=left(-1right)^{n-1}frac{ncdotleft(n+1right)}{2}$

При x=1

$P_n^primeleft(1right)=frac{ncdotleft(n+1right)}{2}$

Для численного интегрирования используются полиномы Лежандра при аргументе принадлежащих отрезку [-1;1]. Это связано с тем, что именно на данном отрезке данные полиномы обладают свойством ортогональности с единичным весом (подробнее смотри мою статью на Дзене «Ортогональные системы функций и ряды Фурье.»). При этом норма квадрата полинома равна $frac{2}{2cdot n+1}$ , что напрямую учитывается при расчете весовых коэффициентов квадратурной формулы.

Численное интегрирование по методу Гаусса‑Лежандра производится по следующей квадратурной формуле:

$int_{a}^{b}{fleft(xright)mathrm{d} x}approxfrac{b-a}{2}cdotsum_{i=0}^{n-1}{w_icdot fleft(frac{b-a}{2}cdot x_i+frac{a+b}{2}right)}$

Где x_i — i‑й корень полинома Лежандра n‑степени; w_i— соответствующие весовые коэффициенты.

Только для полиномов Лежандра 1, 2, 3 и 4 степени можно найти точные значения корней этих полиномов. Для всех полиномов, степени которых больше 4 корни можно найти только численными методами. Для их нахождения рекомендуется использовать рекуррентный метод нахождения корней нелинейных уравнений Исаака Ньютона. (Подробнее смотрите мою статью на Дзене «Решение нелинейных уравнений f(x)=0 методом Ньютона (касательных)»). При начальном приближении, вычисленной по следующей формуле:

$x_i^{left(0right)}=cos{left(picdotfrac{4cdotleft(n-iright)-1}{4cdot n+2}right)}$

1–3 итерации метода Ньютона гарантирует абсолютную точность порядка 10^-10.

Весовые коэффициенты вычисляются по следующей формуле:

$w_i=frac{2}{left(1-x_i^2right)cdotleft(P_n^primeleft(x_iright)right)^2}$

где x_i — i‑й корень полиномов Лежандра n‑й степени.

Программная реализация

По данной методики мною было разработаны на языке C# классы и методы, с помощью которых реализован метод Гаусса‑Лежандра численного интегрирования. Данная реализация позволят вычислять интегралы с использованием данного метода с любым количеством промежуточных точек так как все вычисляется на лету.

Нам понадобится вспомогательный метод, реализованный в классе Equation который ищет приближенно корень нелинейного уравнения методом Ньютона, если известно само уравнение и его производная:

Скрытый текст

        /// <summary>
        /// Вспомогательный метод для оптимального вычисления
        /// i-го корня полиномов Лежандра n-й степени
        /// и соответственного веса для 
        /// интегрирования методом Гаусса-Лежандра 
        /// </summary>
        /// <param name="n">Степень полинома Лежандра</param>
        /// <param name="i">Порядковый номер корня и веса</param>
        /// <returns></returns>
        internal static NodeWeight GetNodeAndWeight(int n, int i)
        {
            double x = Math.Cos(Math.PI * (4.0 * (n - i) - 1.0) / (4.0 * n + 2.0));

            // Встроенный метод Ньютона, оптимизированный под Лежандра
            double x0;
            double p = 0, dp = 0;
            do
            {
                x0 = x;
                // Считаем полином и производную одновременно, за один проход
                p = Polynom(x0, n);
                dp = DeffPolynom(x0, n);
                x = x0 - p / dp;
            }
            while (Math.Abs(x - x0) > 1e-12); // Прецизионная точность

            return new NodeWeight
            {
                Value = x,
                Weight = 2.0 / ((1.0 - x * x) * dp * dp) // Вес считается мгновенно из готовой производной!
            };
        }

Класс LegendrePolynomials реализует все необходимые методы для работы с полиномами Лежандра любой степени, начиная с 3-й:

Скрытый текст

    /// <summary>
    /// Вспомогательная структура для хранения
    /// пары корень полинома Лежандра 
    /// и соответсвующего ему весового коэффициента
    /// </summary>
    internal struct NodeWeight
    {
        internal double Value;  // Корень
        internal double Weight; // Вес
    }
    /// <summary>
    /// Класс для работы с полиномами Лежандра
    /// </summary>
    public static  class LegendrePolynomials
    {
        /// <summary>
        /// Вычисление полинома Лежандра степени n, от аргумента x
        /// </summary>
        /// <param name="x"> аргумент полинома Лежандра </param>
        /// <param name="n"> степень полинома Лежандра</param>
        /// <returns>Вычисленное значение полинома Лежандра</returns>
        public static double Polynom(double x, int n)
        {
            switch(n)
            {
                case 0:
                    {
                        return 1;
                    }
                case 1:
                    {
                        return x;
                    }
                default:
                    {
                        double y1 = 1;
                        double y2 = x;
                        for(int i=2;i<=n;i++)
                        {
                            double y3 = ((2 * i - 1) * x * y2 - (i - 1) * y1) / i;
                            y1 = y2;
                            y2 = y3;
                        }
                        return y2;
                    }
            }
        }
        /// <summary>
        /// Вычисление производной от полинома Лежандра степени n и аргумента x
        /// </summary>
        /// <param name="x">аргумент</param>
        /// <param name="n">степень полинома Лежандра</param>
        /// <returns></returns>
        public static double DeffPolynom(double x, int n)
        {
            // Защита от деления на ноль на краях интервала [-1; 1]
            if (Math.Abs(1.0 - x * x) < 1E-15)
            {
                return (x == 1.0 ? 1.0 : (n % 2 == 0 ? -1.0 : 1.0)) * n * (n + 1) / 2.0;
            }
            return n * (Polynom(x, n - 1) - x * Polynom(x, n)) / (1.0 - x * x);
        }
        /// <summary>
        /// Вычисление i-го корня полинома Лежандра степени n
        /// </summary>
        /// <param name="n">степень полинома </param>
        /// <param name="i"></param>
        /// <returns></returns>
        public static double Root(int n, int i)
        {
            if(n<=i || i<0)
            {
                return double.NaN;
            }
            double x = Math.Cos(Math.PI * (4*(n-i)-1) / (4 * n + 2));
            Func<double, double> f = (double targetX) => Polynom(targetX, n);
            Func<double, double> df = (double targetX) => DeffPolynom(targetX, n);
            return Equation.Nuton(f, df, x, 1e-10);
        }
        /// <summary>
        /// Вычисляется i корень полинома Лежандра n-й степени
        /// </summary>
        /// <param name="n">Степень полинома Лежандра</param>
        /// <param name="i">Порядковый номер корня полинома, который требуется найти</param>
        /// <returns></returns>
        public static double Weight(int n, int i)
        {
            if (n <= i || i < 0)
            {
                return double.NaN;
            }
            double x = Root(n, i);
            double a = DeffPolynom(x, n);
            return 2 / ((1 - x * x) * a * a);
        }
        public static double Weight(int n, double x)
        {
            double a = DeffPolynom(x, n);
            return 2 / ((1 - x * x) * a * a);
        }
        
        /// <summary>
        /// Вспомогательный метод для оптимального вычисления
        /// i-го корня полиномов Лежандра n-й степени
        /// и соответственного веса для 
        /// интегрирования методом Гаусса-Лежандра 
        /// </summary>
        /// <param name="n">Степень полинома Лежандра</param>
        /// <param name="i">Порядковый номер корня и веса</param>
        /// <returns></returns>
        internal static NodeWeight GetNodeAndWeight(int n, int i)
        {
            double x = Math.Cos(Math.PI * (4.0 * (n - i) - 1.0) / (4.0 * n + 2.0));

            // Встроенный метод Ньютона, оптимизированный под Лежандра
            double x0;
            double p = 0, dp = 0;
            do
            {
                x0 = x;
                // Считаем полином и производную одновременно, за один проход
                p = Polynom(x0, n);
                dp = DeffPolynom(x0, n);
                x = x0 - p / dp;
            }
            while (Math.Abs(x - x0) > 1e-12); // Прецизионная точность

            return new NodeWeight
            {
                Value = x,
                Weight = 2.0 / ((1.0 - x * x) * dp * dp) // Вес считается мгновенно из готовой производной!
            };
        }
    }

В классе Integral реализован метод LegendreGause, реализующий численное интегрирование методом Гаусса‑Лежандра.

Скрытый текст

        public static double LegendreGause(Func<double,double> f,double Start, double End, int n=32)
        {
            double Result = 0;
            for(int i=0;i<n;i++)
            {
                NodeWeight node = LegendrePolynomials.GetNodeAndWeight(n,i);
                
                Result += node.Weight * f(((End - Start) * node.Value + End + Start) / 2);
            }
            return (End - Start) * Result / 2;
        }

Тестирование

Тестирование предложенной реализации численного интегрирования будем делать на следующем определенном интеграле:

$int_{0}^{pi}{frac{sin{x}}{1+x}dx}approx0.8438108128003$

Здесь 12 знаков после запятой верны.

Ниже в таблице приводится вычисления данного интеграла по методу Гаусса‑Лежандра и для сравнения по методу Симпсона при разных количествах точках:

Кол‑во точек	Интегрирование методом Гаусса‑Лежандра	Интегрирование методом Симпсона
8	0.8438108709747	0.789880436120
16	0.8438108127999	0.830828283985
32	0.8438108128003	0.840635210390
100		0.843490671272
1000		0.843807632997

Из таблицы видно, что интегрирование по методу Гаусса‑Лежандра уже при 32 точек точность достигается до 12 знаком после запятой, так как метод интегрирования Симпсона при 1000 точках дает всего 5 верных знаков. Это говорит, что метод Гаусса‑Лежандра более эффективен.

Автор: Melamed69

Источник

Информация

Комментарии

Рекомендуем

Динамический метод численного интегрирования Гаусса‑Лежандра на C#: вычисляем узлы и веса «на лету» без переполнения

Теоретическая часть

Программная реализация

Тестирование