Хабр·Russia·Saturday, October 31, 2026 at 12:15 PM UTC

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

Translate

Melamed6910 минут назад

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

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

$P_n\left(x\right)=\frac{1}{2^n}\cdot\sum_{k=0}^{\left[n/2\right]}{\left(-1\right)^k\cdot C_k^n\cdot C_n^{2\cdot n-2\cdot k}\cdot x^{n-2\cdot k}}\,\,\,\,\, (1)$

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

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

$P_0\left(x\right)=1$ $P_1\left(x\right)=x$ $P_{n+1}\left(x\right)=\frac{2\cdot n+1}{n+1}\cdot x\cdot P_n\left(x\right)-\frac{n}{n+1}\cdot P_{n-1}\left(x\right)$

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

$P_n^\prime\left(x\right)=\frac{n}{1-x^2}\cdot\left(P_{n-1}-x\cdot P_n\left(x\right)\right)$

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

При x=-1:

$P_n^\prime\left(-1\right)=\left(-1\right)^{n-1}\frac{n\cdot\left(n+1\right)}{2}$

При x=1

$P_n^\prime\left(1\right)=\frac{n\cdot\left(n+1\right)}{2}$

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

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

$\int_{a}^{b}{f\left(x\right)\mathrm{d} x}\approx\frac{b-a}{2}\cdot\sum_{i=0}^{n-1}{w_i\cdot f\left(\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(0\right)}=cos{\left(pi\cdot\frac{4\cdot\left(n-i\right)-1}{4\cdot n+2}\right)}$

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

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

$w_i=\frac{2}{\left(1-x_i^2\right)\cdot\left(P_n^\prime\left(x_i\right)\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 верных знаков. Это говорит, что метод Гаусса-Лежандра более эффективен.

Если эта публикация вас вдохновила и вы хотите поддержать автора — не стесняйтесь нажать на кнопку

View the original on Хабр →

KioskNews shows a cleaned-up reading view extracted from the publisher’s page — the original always lives on their site, not ours.