CC BY
216
УДК 519.622.2
Д.М. Доронин
ЧИСЛЕННЫЙ МЕТОД ИНТЕГРИРОВАНИЯ НА ОСНОВЕ АППРОКСИМАЦИИ ЛАГРАНЖА И КВАДРАТУРНЫХ ФОРМУЛ ГАУССА-КРИСТОФФЕЛЯ
Описан алгоритм численного метода решения обыкновенных дифференциальных уравнений. Метод позволяет достаточно просто получать как явные, так и неявные многошаговые алгоритмы на основе аппроксимации Лагранжа и квадратурных формул Гаусса-Кристоффеля. Показано, что явный метод второго, третьего и четвёртого порядков совпадает с соответствующими формами ме-
тода Адамса-Башфорта. Предложена новая численная схема интегрирования высокого порядка, улучшающая сходимость описанного метода.
Задача Коши, дифференциальные уравнения, аппроксимация Лагранжа, формулы Г аусса-Кристоффеля
D.M. Doronin A METHOD FOR NUMERICAL INTEGRATION BASED ON LAGRANGE APPROXIMATION AND GAUSS-CHRISTOFFEL QUADRATURE FORMULAS
The article presents a numerical method offered for the ordinary differential equations. The method allows getting explicit and implicit multistep dynamic-order algorithms on the base of Lagrange Approximation and Gauss-Christoffel quadrature formulas. It is founded that the second, the third, and the forth order explicit methods coincide with the Adams-Bashforth equations. A new numerical schema was proposed to improve the convergence of the described method.
A Cauchy problem, differential equations, Lagrange polynomials, Gauss-Christoffel equations
Введение
В научных приложениях важной задачей является улучшение сходимости разностных схем, используемых при решении систем дифференциальных уравнений [1-3]. Однако более высокий порядок схемы приводит к большему числу операций, связанных с вычислением подынтегральной функции, что в свою очередь приводит к существенному снижению быстродействия программы. Таким образом, актуальной является задача разработки эффективного метода интегрирования, который можно использовать при решении систем обыкновенных дифференциальных уравнений.
В работе предлагается аппроксимировать подынтегральную функцию полиномом Лагранжа и использовать формулы наивысшей алгебраической точности (Гаусса-Кристоффеля) [4] для вычисления значения интеграла. Следует заметить, что алгоритм получения численных схем интегрирования принципиально отличается от классических методов Эйлера, Рунге-Кутты и Адамса, поскольку в методах Эйлера и Рунге-Кутты используется разложение в ряд Тейлора, а в методе Адамса - полином Ньютона.
Постановка задачи
Пусть имеется обыкновенное дифференциальное уравнение N-го порядка (1). Это уравнение, как указано в [1], с помощью замены (2) можно свести к эквивалентной системе N дифференциальных уравнений первого порядка (3).
u(N) (x) = f (x, u, u ', u ',..., u( N-1)) (1)
uk (x) = u(k)(x); u'k (x) = uk+1(x); 1 < k < N (2)
u ' (x) = u{;
u" (x )= u 2; (3)
u( N) (x ) = f ( x, u, u ' , u '”,..., u( N-1)).
Необходимо найти функцию u(x) при заданных начальных условиях (4). Это так называемая задача Коши [2].
u(x0 ) = u(0);
u' (x0 ) = u(1);
< u'1 (x0 ) = u(2); (4)
u'N-1 (x0 ) = u (N}.
Нас будут интересовать системы обыкновенных дифференциальных уравнений, для которых выполняются условия теоремы Пикара [2] и, следовательно, решение задачи Коши существует и единственно.
12
Описание метода решения системы дифференциальных уравнений на основе аппроксимации Лагранжа
В работе предложен метод решения системы уравнений (3), в соответствии с которым подынтегральная функция аппроксимируется полиномом в форме Лагранжа. Далее для вычисления полученного интеграла используются формулы Гаусса-Кристоффеля, известные как формулы наивысшей алгебраической точности.
Будем решать систему уравнений (3) предложенным методом. Для этого рассмотрим задачу Коши для одного уравнения первого порядка (5), затем обобщим решение на случай системы.
и '(*) = I ( х, и( х)), (5)
где переменная х определена на некотором отрезке [а; Ь].
Интегрируя уравнение (5), преобразуем его к эквивалентному интегральному уравнению типа Вольтерра (6).
X
и(х) = и(а) +1 / (т, и(х))Ст (6)
а
Обозначим точное решение уравнения переменной у и разобьём исследуемую область [а; Ь] на N частей. Обозначим значение аргумента х на левой границе х0=а, на правой границе - х^=Ь, тогда для произвольного узла с индексом i е [0; N] будет справедливо выражение
х = хо + (Ь - а) N (7)
Далее будем рассматривать подынтегральную функцию в (6) не на всей плоскости её аргументов, а на определённой интегральной кривой и(х), соответствующей искомому решению, тогда она будет функцией одного аргумента
Р (х) = I (х, и ( х)) (8)
Пусть известно приближённое решение уп в некоторой точке хп, тогда приближённое решение в следующей точке сетки будет определяться по формуле
хп+1
Уп+1 (х)= Уп + | Р(х)Сх (9)
х
Аппроксимируем подынтегральную функцию полиномом Лагранжа к-го порядка по к известным приближённым решениям в предыдущих точках. Тогда уравнение (9) будет иметь следующий вид:
п хп+1
Уп+1 = Уп + X Р (х* ) 11 (х)Сх , ( 10)
п
г=п-к
где I (х) = П — —— - базис полинома Лагранжа.
к=0,jфi х} — хк
Далее линейным преобразованием аргумента (11) перейдём к пределам интегрирования [-1; 1], при которых подынтегральная функция представляет собой формулу Гаусса-Кристоффеля с весовым множителем, равным единице. Значение подынтегральной функции находится с помощью узлов и весов полинома Лежандра [5], обозначенных переменными Хк и Ук соответственно. Обратным линейным преобразованием можно получить узлы и веса для произвольного отрезка [а, Ь].
хк = 2(а+Ь)+2(Ь—а)Х; ск=\(Ь—а)?к (11) Узлы полиномов Лежандра Хк могут быть определены, например, одним из методов одномерной оптимизации, а веса ук по формуле (12), результаты приведены в табл. 1.
*=ик! СЬ' (|2)
1 сп (х2 — 1)п
где Р (х) =----------------полином Лежандра степени п.
п 2пп! Схп
х
Таблица 1
Узлы Х и веса ук полиномов Лежандра. к - индекс коэффициента, М - порядок полинома Лежандра
Используя формулы Гаусса-Кристоффеля, получим выражение (13) для приближённого вычисления подынтегральной функции в (10).
Ь М
11,(х)Сх @ Xcjl^ (xj), (13)
а
где с . = — (Ь — а)Г., М - порядок полинома Лежандра, х . = —(а + Ь) +—(Ь — а)£..
; 2 ; ; 2 2 ;
Для равномерного шага а = хп, Ь = хп+1, Ь — а = к, а + Ь = 2хп + к с учётом (14) получим
формулу (15).
к к .
с. = 2 ; х. = хп + 2(1 + Х) (14)
И п м
Уп+1 = Уп + - Ёр(-)ЁУ(-«У), (15)
2 г=п-к у=1
X 1
где О у = +—(1 + ^; ) .
И 2
Запишем базисную функцию Лагранжа в более удобной для численных вычислений форме (16).
п X - X п ^ + (П - У)
1 (-)= П — = П — (—б)
у=п-к, - - -у ]=п-к, 1 - У
Точность численных вычислений (16) может быть повышена, если использовать не вещественные значения -к, а их целочисленные индексы к
. . п -(1 + Ху) + (п - У)
I,(-«, )= П --—- (17)
У=п-к, 1 у
1
Поскольку в выражение (17) входят константные значения, не зависящие ни от пределов интегрирования, ни от текущего значения переменного аргумента, его можно преобразовать к виду (18). Значения коэффициентов Л,j для разных порядков полинома Лагранжа представлены в табл. 2, для вычисления последних использовались полиномы Лежандра 5-го порядка:
И п
Уп+1 = Уп + — (- ) (18)
24 г=п-к
Таблица 2
Значения коэффициентов Л, для разных порядков полинома Лагранжа
к Лп 8 п-8 Лп-7 Лп-6 Лп-5 Лп-4 Лп-3 Лп-2 Лп-! Лп
1 0 0 0 0 0 0 0 36 -12
2 0 0 0 0 0 0 46 -32 10
3 0 0 0 0 0 55 -59 37 -9
4 0 0 0 0 63 -92 87 -42 8.37
5 0 0 0 71 -132 166 -122 48 -7.92
6 0 0 79 -177 280 -273 162 -53 7.57
7 0 86 -229 433 -529 417 -207 59 -7.3
8 93 -285 631 -925 912 -603 257 -64 7.08
Было замечено, что полученные значения Лi совпадают с коэффициентами метода Адамса
соответствующего порядка решения систем обыкновенных дифференциальных уравнений. Представленный в работе метод может быть полезен для получения коэффициентов численных схем Адамса-
Башфорта и Адамса-Башфорта-Мултона высокого порядка, аналитическое вычисление которых вызывает определённые трудности.
Оценка эффективности метода для жёстких задач
Оценка эффективности предложенного метода осуществлялась на модельной задаче Коши (19) с отрицательными коэффициентами ai, определяющими жёсткость задачи.
^ = агу + £ (-), 1 <, < п, 0 < - < 1, (19)
дх
Для данной задачи можно ввести коэффициент жёсткости (20).
шах|Яе(- Л, (-))
5 (х) = ш<п' -^ (20)
шт Яе(- Л, (х))
1<1<п' '
где Лi - собственные числа якобиана системы уравнений (19).
Систему дифференциальных уравнений считают жёсткой, если коэффициент жёсткости больше единицы 5(х) > 1 .
Оценим эффективность применения разработанного метода различных порядков на примере задачи (21) [5], являющейся частным случаем системы уравнений (19).
Фі =
йх
йу2 = йх
998уі + 1998у 29 уі(0) = 1
-999уі -1999у2,у2(0) = 0
(21)
Собственные значения якобиана 3 таковы: Л1 = 2560 и Л2 = 1559 - следовательно, коэффициент жёсткости 5 = 1.64 . Согласно [5] система является жёсткой, а её точное решение (21) может быть найдено по аналитической формуле (22).
(22)
\ У1(х) = 2ехр(-4х) - ехр(-10000х) [ у 2( х) = - ехр(-4 х) - ехр(-10000х)
На рисунке приведена оценка отклонения численного решения системы уравнений (21) от точного решения (22) на отрезке 0 £ х £ 0.01. Установлено, что оптимальный результат наблюдается для схем 6-го порядка, когда к = 5, далее с увеличением порядка ошибка растёт. Это объясняется ухудшением интерполяции на границах интервала интегрирования. Для устранения этого дефекта предлагается выбирать узлы интерполяции в соответствии с нулями полинома Чебышева. Такая методика приводит к новым схемам интегрирования, в частности для метода восьмого порядка была получена схема
к
у„, = у, ^(52.298^,-6 -41.820^,-5 + 21.2/ -12.420/,-, + 4.698/,)
24
(23)
Распределение отклонения численного решения системы дифференциальных уравнений от точного решения для разработанного метода различных порядков (для К=7 представлена модифицированная схема)
Заключение
Получен новый многошаговый численный метод интегрирования на основе аппроксимации Лагранжа и квадратурных формул Гаусса-Кристоффеля для решения задач Коши и численного вычисления интегралов. Метод позволяет достаточно просто получать явные и неявные многошаговые алгоритмы без обычно используемых громоздких выводов. Показано, что явные методы второго, третьего и четвёртого порядков совпадают с соответствующими формами метода Адамса-Башфорта, вероятно, разработанный метод будет совпадать со схемами Адамса и для более высоких порядков. Проведённый анализ устойчивости метода при решении жёстких задач показал, что явные схемы метода не являются достаточно устойчивыми, поэтому для этого класса задач интерес представляют неявные схемы. Предложена методика улучшения качества схем высокого порядка за счёт выбора
узлов интерполяции вблизи нулей полиномов Чебышева, представлена модифицированная схема явного метода восьмого порядка. Для решения прикладных задач рекомендуется использовать явную схему шестого порядка или модифицированную схему восьмого порядка.
ЛИТЕРАТУРА
1. Каханер Д. Численные методы и программное обеспечение / Д. Каханер, К. Моулер, С. Нэш. М.: Мир, 2001. 575 с.
2. Калиткин Н.Н. Численные методы / Н.Н. Калиткин. М.: BHV, 2011. 596 с.
3. Хайрер Э. Решение обыкновенных дифференциальных уравнений. Жёсткие и дифференциально-алгебраические задачи / Э. Хайрер, Г. Ваннер. М.: Мир, 1999. 685 с.
4. Самарский А.А. Численные методы / А.А. Самарский, А.В. Гулин. М.: Наука, 1989. 432 с.
5. Numerical Recipes in C : the art of scientific computing / William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery. 2nd ed. Cambridge University Press, 1993. 1020 p.
Доронин Дмитрий Михайлович - Dmitri M. Doronin -
аспирант кафедры «Радиотехника graduate department of Radiotrician
и электродинамика» Саратовского and electrodynamics
государственного университета N.G. Chernyshevski Saratov State University
имени Н.Г. Чернышевского
Статья поступила в редакцию 10.11.13, принята к опубликованию 15.12.13
