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

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

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

Полиномы Лежандра 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 верных знаков. Это говорит, что метод Гаусса-Лежандра более эффективен.

ссылка на оригинал статьи https://habr.com/ru/articles/1102060/