close

Вход

Забыли?

вход по аккаунту

?

Численный метод интегрирования на основе аппроксимации Лагранжа и квадратурных формул Гаусса-Кристоффеля.

код для вставкиСкачать
Математика и механика
УДК 519.622.2
Д.М. Доронин
ЧИСЛЕННЫЙ МЕТОД ИНТЕГРИРОВАНИЯ НА ОСНОВЕ АППРОКСИМАЦИИ
ЛАГРАНЖА И КВАДРАТУРНЫХ ФОРМУЛ ГАУССА-КРИСТОФФЕЛЯ
Описан алгоритм численного метода решения обыкновенных дифференциальных уравнений. Метод позволяет достаточно просто получать как явные, так
и неявные многошаговые алгоритмы на основе аппроксимации Лагранжа и квадратурных формул Гаусса-Кристоффеля. Показано, что явный метод второго,
третьего и четвёртого порядков совпадает с соответствующими формами ме11
Вестник СГТУ. 2013. № 4 (73)
тода Адамса-Башфорта. Предложена новая численная схема интегрирования высокого порядка, улучшающая сходимость описанного метода.
Задача Коши, дифференциальные уравнения, аппроксимация Лагранжа,
формулы Гаусса-Кристоффеля
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, GaussChristoffel equations
Введение
В научных приложениях важной задачей является улучшение сходимости разностных схем,
используемых при решении систем дифференциальных уравнений [1-3]. Однако более высокий порядок схемы приводит к большему числу операций, связанных с вычислением подынтегральной
функции, что в свою очередь приводит к существенному снижению быстродействия программы. Таким образом, актуальной является задача разработки эффективного метода интегрирования, который
можно использовать при решении систем обыкновенных дифференциальных уравнений.
В работе предлагается аппроксимировать подынтегральную функцию полиномом Лагранжа и использовать формулы наивысшей алгебраической точности (Гаусса-Кристоффеля) [4] для вычисления значения интеграла. Следует заметить, что алгоритм получения численных схем интегрирования принципиально отличается от классических методов Эйлера, Рунге-Кутты и Адамса, поскольку в методах Эйлера и
Рунге-Кутты используется разложение в ряд Тейлора, а в методе Адамса – полином Ньютона.
Постановка задачи
Пусть имеется обыкновенное дифференциальное уравнение N-го порядка (1). Это уравнение,
как указано в [1], с помощью замены (2) можно свести к эквивалентной системе N дифференциальных уравнений первого порядка (3).
u ( N ) (x ) = f ( x, u, u ′, u ′′,..., u ( N −1) )
(1)
(k )
(2)
u k ( x ) = u ( x ); u k′ (x ) = u k +1 ( x ); 1 ≤ k ≤ N
u ' ( x ) = u1 ;
u ' ' (x ) = u ;

2
(3)

...

u ( N ) (x ) = f ( x, u , u ′, u ′′,..., u ( N −1) ).
Необходимо найти функцию u(x) при заданных начальных условиях (4). Это так называемая
задача Коши [2].
u ( x0 ) = u ( 0) ;

(1)
u ' ( x0 ) = u ;

( 2)
(4)
u '1 ( x0 ) = u ;
...

u ' N −1 ( x0 ) = u ( N ) .

Нас будут интересовать системы обыкновенных дифференциальных уравнений, для которых
выполняются условия теоремы Пикара [2] и, следовательно, решение задачи Коши существует и
единственно.
12
Математика и механика
Описание метода решения системы дифференциальных уравнений на основе аппроксимации Лагранжа
В работе предложен метод решения системы уравнений (3), в соответствии с которым подынтегральная функция аппроксимируется полиномом в форме Лагранжа. Далее для вычисления полученного интеграла используются формулы Гаусса-Кристоффеля, известные как формулы наивысшей
алгебраической точности.
Будем решать систему уравнений (3) предложенным методом. Для этого рассмотрим задачу
Коши для одного уравнения первого порядка (5), затем обобщим решение на случай системы.
u ′( x ) = f ( x, u ( x)) ,
(5)
где переменная x определена на некотором отрезке [a; b].
Интегрируя уравнение (5), преобразуем его к эквивалентному интегральному уравнению типа
Вольтерра (6).
x
u ( x ) = u (a) + ∫ f (τ , u ( x))dτ
(6)
a
Обозначим точное решение уравнения переменной y и разобьём исследуемую область [a; b]
на N частей. Обозначим значение аргумента x на левой границе x0=a, на правой границе – xN=b, тогда
для произвольного узла с индексом i ∈ [0; N ] будет справедливо выражение
xi = x0 + (b − a )
i
N
(7)
Далее будем рассматривать подынтегральную функцию в (6) не на всей плоскости её аргументов, а на определённой интегральной кривой u(x), соответствующей искомому решению, тогда
она будет функцией одного аргумента
(8)
F ( x ) = f ( x, u ( x ))
Пусть известно приближённое решение yn в некоторой точке xn, тогда приближённое решение
в следующей точке сетки будет определяться по формуле
y n +1 (x ) = y n +
x n +1
∫ F ( x)dx
(9)
xn
Аппроксимируем подынтегральную функцию полиномом Лагранжа k-го порядка по k известным
приближённым решениям в предыдущих точках. Тогда уравнение (9) будет иметь следующий вид:
y n +1 = y n +
где l j ( x) =
n
∏
k =0 , j ≠i
xn +1
n
∑ F ( x ) ∫ l ( x)dx ,
i=n−k
i
i
(10)
xn
x − xk
– базис полинома Лагранжа.
x j − xk
Далее линейным преобразованием аргумента (11) перейдём к пределам интегрирования [-1; 1], при
которых подынтегральная функция представляет собой формулу Гаусса-Кристоффеля с весовым
множителем, равным единице. Значение подынтегральной функции находится с помощью узлов и
весов полинома Лежандра [5], обозначенных переменными ξ k и γ k соответственно. Обратным линейным преобразованием можно получить узлы и веса для произвольного отрезка [a, b].
1
1
1
(11)
xk = (a + b) + (b − a )ξ k ; ck = (b − a )γ k
2
2
2
Узлы полиномов Лежандра ξ k могут быть определены, например, одним из методов одномерной оптимизации, а веса γ k по формуле (12), результаты приведены в табл. 1.
PM ( χ )
1
γk =
dχ ,
∫
PM′ (ξ k ) −1 χ − ξ k
1
(12)
1 d n ( x 2 − 1) n
– полином Лежандра степени n.
где Pn ( x) = n
2 n!
dx n
13
Вестник СГТУ. 2013. № 4 (73)
Узлы
k
3
4
5
6
7
8
9
и веса
Таблица 1
γk
– индекс коэффициента,
полиномов Лежандра.
M
– порядок полинома Лежандра
k
1
2
ξk
-0.5774
0.5774
γk
1.0000
1.0000
ξk
-0.7746
0.0000
0.7746
γk
0.5556
0.8889
0.5556
ξk
-0.8611
-0.3400
0.3400
0.8611
γk
0.3479
0.6521
0.6521
0.3479
ξk
-0.9062
-0.5385
0.0000
0.5385
0.9062
γk
0.2369
0.4786
0.5689
0.4786
0.2369
ξk
-0.9325
-0.6612
-0.2386
0.2386
0.6612
0.9325
γk
0.1713
0.3608
0.4679
0.4679
0.3608
0.1713
ξk
-0.9491
-0.7415
-0.4058
0.0000
0.4058
0.7415
0.9491
γk
0.1295
0.2797
0.3818
0.4180
0.3818
0.2797
0.1295
ξk
-0.9603
-0.7967
-0.5255
-0.1834
0.1834
0.5255
0.7967
0.9603
γk
0.1012
0.2224
0.3137
0.3627
0.3627
0.3137
0.2224
0.1012
ξk
-0.9682
-0.8360
-0.6134
-0.3243
0.0000
0.3243
0.6134
0.8360
0.9682
γk
0.0813
0.1806
0.2606
0.3123
0.3302
0.3123
0.2606
0.1806
0.0813
M
2
ξk
3
4
5
6
7
8
9
Используя формулы Гаусса-Кристоффеля, получим выражение (13) для приближённого вычисления подынтегральной функции в (10).
b
M
a
j =1
∫ li ( x)dx ≅ ∑ c j li ( x j ) ,
(13)
1
1
1
(b − a)γ j , M – порядок полинома Лежандра, x j = (a + b) + (b − a)ξ j .
2
2
2
Для равномерного шага a = xn , b = xn+1 , b − a = h , a + b = 2 xn + h с учётом (14) получим
где c j =
формулу (15).
cj =
14
h
h
γ j ; x j = xn + (1 + ξ j )
2
2
(14)
Математика и механика
y n+1 = y n +
где Ω j =
M
h n
F ( xi )∑ γ j li (hΩ j ) ,
∑
2 i = n −k
j =1
(15)
xn 1
+ (1 + ξ j ) .
h 2
Запишем базисную функцию Лагранжа в более удобной для численных вычислений форме (16).
x − xn
+ (n − j )
n
n
x − xj
h
(16)
li ( x ) = ∏
= ∏
i− j
j = n −k , xi − x j
j = n −k ,
i≠ j
i≠ j
Точность численных вычислений (16) может быть повышена, если использовать не вещественные значения xk, а их целочисленные индексы k
1
(1 + ξ j ) + (n − j )
li (hΩ j ) = ∏ 2
i− j
j =n−k ,
n
(17)
i≠ j
Поскольку в выражение (17) входят константные значения, не зависящие ни от пределов интегрирования, ни от текущего значения переменного аргумента, его можно преобразовать к виду (18).
Значения коэффициентов λi для разных порядков полинома Лагранжа представлены в табл. 2, для
вычисления последних использовались полиномы Лежандра 5-го порядка:
yn+1 = yn +
Значения коэффициентов
λi
h n
∑ λi F ( xi )
24 i =n−k
(18)
Таблица 2
для разных порядков полинома Лагранжа
k
λn−8
λ n −7
λ n −6
λn−5
λn − 4
λn−3
λn − 2
λn −1
λn
1
2
3
4
5
6
7
8
0
0
0
0
0
0
0
93
0
0
0
0
0
0
86
-285
0
0
0
0
0
79
-229
631
0
0
0
0
71
-177
433
-925
0
0
0
63
-132
280
-529
912
0
0
55
-92
166
-273
417
-603
0
46
-59
87
-122
162
-207
257
36
-32
37
-42
48
-53
59
-64
-12
10
-9
8.37
-7.92
7.57
-7.3
7.08
Было замечено, что полученные значения λi совпадают с коэффициентами метода Адамса
соответствующего порядка решения систем обыкновенных дифференциальных уравнений. Представленный в работе метод может быть полезен для получения коэффициентов численных схем АдамсаБашфорта и Адамса-Башфорта-Мултона высокого порядка, аналитическое вычисление которых вызывает определённые трудности.
Оценка эффективности метода для жёстких задач
Оценка эффективности предложенного метода осуществлялась на модельной задаче Коши
(19) с отрицательными коэффициентами ai , определяющими жёсткость задачи.
dyi
= ai y + f i ( x), 1 ≤ i ≤ n, 0 < x ≤ 1,
dx
Для данной задачи можно ввести коэффициент жёсткости (20).
max Re (− λi ( x ) )
S ( x ) = 1≤i ≤n
min Re(− λi ( x ) )
(19)
(20)
1≤i ≤ n
где λi – собственные числа якобиана системы уравнений (19).
Систему дифференциальных уравнений считают жёсткой, если коэффициент жёсткости
больше единицы S ( x) > 1 .
Оценим эффективность применения разработанного метода различных порядков на примере
задачи (21) [5], являющейся частным случаем системы уравнений (19).
15
Вестник СГТУ. 2013. № 4 (73)
 dy1
 dx = 998 y1 + 1998 y 2 , y1 (0) = 1
,
(21)

 dy2 = −999 y − 1999 y , y (0) = 0
1
2
2
 dx
Собственные значения якобиана J таковы: λ1 = 2560 и λ2 = 1559 – следовательно, коэффициент жёсткости S = 1.64 . Согласно [5] система является жёсткой, а её точное решение (21) может
быть найдено по аналитической формуле (22).
 y1 ( x) = 2 exp(−4 x ) − exp(−10000 x)
(22)

 y 2 ( x ) = − exp(−4 x ) − exp(−10000 x )
На рисунке приведена оценка отклонения численного решения системы уравнений (21) от
точного решения (22) на отрезке 0 ≤ x ≤ 0.01 . Установлено, что оптимальный результат наблюдается
для схем 6-го порядка, когда k = 5 , далее с увеличением порядка ошибка растёт. Это объясняется
ухудшением интерполяции на границах интервала интегрирования. Для устранения этого дефекта
предлагается выбирать узлы интерполяции в соответствии с нулями полинома Чебышева. Такая методика приводит к новым схемам интегрирования, в частности для метода восьмого порядка была
получена схема
yn+1 = yn +
h
(52.298 f n−6 − 41.820 f n−5 + 21.244 f n−3 − 12.420 f n−1 + 4.698 f n )
24
(23)
Распределение отклонения численного решения системы дифференциальных уравнений от точного решения
для разработанного метода различных порядков (для K=7 представлена модифицированная схема)
Заключение
Получен новый многошаговый численный метод интегрирования на основе аппроксимации
Лагранжа и квадратурных формул Гаусса-Кристоффеля для решения задач Коши и численного вычисления интегралов. Метод позволяет достаточно просто получать явные и неявные многошаговые
алгоритмы без обычно используемых громоздких выводов. Показано, что явные методы второго, третьего и четвёртого порядков совпадают с соответствующими формами метода Адамса-Башфорта, вероятно, разработанный метод будет совпадать со схемами Адамса и для более высоких порядков.
Проведённый анализ устойчивости метода при решении жёстких задач показал, что явные схемы метода не являются достаточно устойчивыми, поэтому для этого класса задач интерес представляют
неявные схемы. Предложена методика улучшения качества схем высокого порядка за счёт выбора
16
Математика и механика
узлов интерполяции вблизи нулей полиномов Чебышева, представлена модифицированная схема явного метода восьмого порядка. Для решения прикладных задач рекомендуется использовать явную
схему шестого порядка или модифицированную схему восьмого порядка.
ЛИТЕРАТУРА
1. Каханер Д. Численные методы и программное обеспечение / Д. Каханер, К. Моулер, С.
Нэш. М.: Мир, 2001. 575 с.
2. Калиткин Н.Н. Численные методы / Н.Н. Калиткин. М.: BHV, 2011. 596 с.
3. Хайрер Э. Решение обыкновенных дифференциальных уравнений. Жёсткие и дифференциально-алгебраические задачи / Э. Хайрер, Г. Ваннер. М.: Мир, 1999. 685 с.
4. Самарский А.А. Численные методы / А.А. Самарский, А.В. Гулин. М.: Наука, 1989. 432 c.
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
17
Документ
Категория
Без категории
Просмотров
12
Размер файла
635 Кб
Теги
численные, лагранжа, метод, формула, основы, кристоффеля, интегрированный, квадратурные, гаусса, аппроксимация
1/--страниц
Пожаловаться на содержимое документа