Slipm - программа на языке Maple для численного решения частичной проблемы Штурма-Лиувилля на основе непрерывного аналога метода Ньютона Текст научной статьи по специальности «Математика»

УДК 517.9, 519.6, 539.12

SLIPM — программа на языке MAPLE для численного решения частичной проблемы Ш^турма—Лиувилля на основе непрерывного аналога метода Ньютона

И. В. Пузынин, Т. П. Пузынина, В. Ч. Тхак

Лаборатория информационных технологий Объединённый институт ядерных исследований, Дубна, Россия Laboratory of Information Technologies Joint Institute for Nuclear Research 141980, Dubna, Russia

SLIPM — (The Sturm-LIouville Problem in Maple) — программа на языке системы компьютерной алгебры (СКА) MAPLE, предназначенная для численного решения с помощью непрерывного аналога метода Ньютона (НАМН) частичной проблемы Штурма— Лиувилля, т.е. для вычисления некоторого собственного значения линейного дифференциального оператора второго порядка и соответствующей собственной функции, удовлетворяющей однородным граничным условиям. SLIPM является развитием написанных на языке Фортран программ SLIP1 и SLIPH4 и дополнена двумя новыми способами вычисления начального значения итерационного параметра то.

Ключевые слова: вычислительная математика, непрерывные аналоги итерационных методов, системы компьютерной алгебры, MAPLE.

1. Постановка задачи

Для линейного дифференциального уравнения

^(1)(А, у) = у" (х) + 2р(х)у'(х) + (q(x) — Хг(х))у(х) = 0, а < х < Ь, (1)

требуется вычислить собственное значение А и соответствующее собственное решение у(х), удовлетворяющее однородным краевым условиям

р(2) (А, у) = di (X, а)у'(а) + fi (X, а)у(а) = 0, (2)

V(3)(X, у) = d2(X, Ъ)у'(Ъ) + f2(X, Ь)у(Ь) = 0. (3)

Для исключения тривиальных решений у(х) = 0 вводится условие нормировки

ь

>/4)(Х,у) = j y2(x)dx = 1. (4)

а

Здесь р(х), q(x), г(х) — заданные функции, обеспечивающие существование нетривиальных собственных решений у = уп(х) граничной задачи (1)—(4), которым соответствуют собственные значения параметра А = Ап. Функции dj, fj непрерывно дифференцируемы по А и d2 + fj > 0; j = 1, 2.

Программа SLIPM по заданному начальному приближению {ХП,уП} вычисляет частичное решение {ХП,уП} задачи (1)-(4), где п — фиксировано, х G ши, Wh — равномерная разностная сетка на отрезке а^ х ^ b с шагом h. Алгоритмы численного решения задачи, взятые нами за основу, реализованы в SLIP1 [1] и SLIPH4 [2] и хорошо зарекомендовали себя при решении многих сложных физических задач [3,4].

Статья поступила в редакцию 28 ноября 2009 г.

Работа выполнена при финансовой поддержке РФФИ, грант 09-01-00770.

2. Метод решения (НАМН)

2.1. Алгоритм вычисления начального приближения

Алгоритм вычисления начального приближения основан на решении двух «встречных» задач Коши. Для первой задачи Коши (х ^ а) ставятся начальные условия в точке х = а, где условие (2) дополняется условием на производную, если = 0. Аналогичное условие для второй задачи (х ^ Ь) ставится в точке х = Ь. Задача об отыскании собственных значений сводится тогда к отысканию корней функции

т(^)= у'п(хт)лев * уп(хт)прав - уп(хт)лев * у'п(хт) прав = 0, (5)

появляющейся из условия равенства логарифмических производных для решений двух задач Коши в некоторой внутренней точке отрезка хт € [а,Ь]

у'п(хт)лев _ Уп(хт)прав Уп(хт)лев уп(хт)прав'

(6)

Если задачи Коши решать приближённо с помощью трёхточечной разностной схемы точности 0(Ь2), где к — шаг равномерной сетки шь, а функции dj(X), ¡^(X) являются полиномами от степени Л, что справедливо для широкого круга практических задач, то функция Т(X) является полиномом. Его корни Хп находятся по методу Ньютона с исключением уже найденных корней [5,6] (к — номер итерации):

\&+1 = \к__Т (^п)__(7)

" = " Ц№) - Т(Л^)Е^-11 Ф;' ()

Итерационный процесс (7) прекращается при одновременном выполнении двух условий:

|А^+1 - Хкп1 < ет и |Г(Хкп+1)1 < ет, (8)

где ет > 0 — заданное малое число. После определения \п с заданной точностью мы одновременно получаем и приближенное значение соответствующей сеточной собственной функции уп(х^). В БЫРМ этот алгоритм реализован в процедуре Ж^^ГОШ.

2.2. Итерационные схемы на основе НАМН

Задача (1)—(4) является нелинейным функциональным уравнением

р(г) = 0, г = (X,y(x)),

где (р — совокупность функций , ] = 1,2, 3,4. Согласно НАМН [3,4], задача (1)—(4) заменяется эволюционным уравнением

у'(г(г))г'(г) = -у(г(г)), 0 < ж, (9)

с начальным условием

*(0) = ^ = (Хо,уо(х)). (10)

Здесь — производная Фреше оператора <р, г'(£) = (у,(£),у(х,£)), где

М*) = \'(1), у(х,1) = у' (х,1). (11)

При достаточно общих предположениях на р(£) имеем — 2* \\ = 0, где

= (\*,у*(х)) — решение задачи (1)-(4), ближайшее к г0. Для приближенного решения задачи (9)-(10) методом Эйлера введём дискретное множество узлов

{^к} = [тк = tk+l — Ьк ,к = 0,1,2,...} и запишем приближенные выражения для соотношений (11)

' ау(х,ь)

АЩ)

)

= у{х,ьк+1) — у(х,гк)

Тк ''

х^к+1) — )

Тк

(12)

Введём обозначения Хк = Х^к), Ук = у(х,Ьк), Цк = ^(Ьк), Vк = у(х,Ьк). При дискретной аппроксимации по переменной х эволюционного уравнения (9) использовались трёхточечные разностные и квадратурные формулы с точностью аппроксимации порядка 0(Ь?) на равномерной сетке шь € [а,Ь]. Для вычисления Хк+1,Ук+1 при известных Хк,ук на каждом шаге итерации необходимо: 1) решить относительно тк уравнение

(Хк ,Ук)шк = —фх'ук, 3 = 1, 2, 3; 2) вычислить значение Цк из уравнения ( (■, ■) = ■ )ёж)

^к

1 + (Ук,Ук),

2(ук ,шк)

3) найти новые приближения А&+1 и Ук+1, используя (12), т.е.

Хк+1 = Хк + Тк^к, Ук+1 = (1 — тк )ук + тк цк Шк,

(13)

(14)

(15)

где Тк определяется специальными алгоритмами. Итерации прекращаются при выполнении условия

4 < £, (16) где 5к определяется по одной из формул:

5к = тахтах 1^(з'(Хк,Ук(ж))|, 3 = 1, 2, 3,

з х

(17)

5к =

<р(1' (Хк,ук(х))Ах

(18)

х € шь, £ > 0 — заданное малое число.

2.3. Алгоритмы вычисления итерационного параметра т^

Вычисление итерационного параметра Тк связано с изменением невязки 5к в ходе итераций. Ниже приводятся следующие алгоритмы выбора Тк, в которых т0 — некоторое заданное значение 0 <т0 ^ 1.

1. Тк = то.

(19)

Этот алгоритм при достаточно малом то обычно применяется при плохих начальных приближениях с целью проверить возможность сходимости от этих приближений. Сходимость при этом очень медленная. При Тк = 1 получается классическая схема Ньютона.

2. Тк

тт(1,2тк-1), если 5к < 5к-1, тах(то,Тк-1/2), если 5к > 5к-1,

ь

где 5к определяется по формуле (17).

Этот алгоритм аналогичен широко распространённому способу выбора шага интегрирования в стандартных программах решения задачи Коши, вычисления интегралов.

Алгоритм рекомендуется применять при хороших начальных приближениях. Он обеспечивает быструю сходимость, однако не всегда устойчив.

-к-1 ' - -к '

тах(то,тк-1 ), если 5к > бк-\,

, тт(1,тк-1 Л-1), если 5к < 5к-\, .

Ъ.^ =\ _______ -к-1\ - - с . г (21)

где 5к определяется по формуле (17). Этот алгоритм более устойчив и обеспечивает сходимость в достаточно широкой области начальных приближений, однако зависит от задания т0.

4 Тк = (22)

ёк-1 + 5к (1)'

где 5к определяется по формуле (18), 5к(1) — невязка на к— той итерации для Тк = 1. Это алгоритм оптимального выбора Тк [7].

5. На равномерной сетке шт на отрезке [0,1] с шагом Ат вычисляется такое значение Тк, которому соответствует минимальная невязка. Этот алгоритм более общий, чем в п.4, хотя и требует большего объёма вычислений.

Недостатком алгоритмов 1-3 является эмпирическое задание начального значения параметра т0.

2.4. Алгоритмы вычисления значения параметра т0

В данной работе предложены и программно реализованы два новых алгоритма вычисления 70:

1 . - = 2^), (23)

II. то = 2

So

(24)

где невязка 5к определяется по формуле (17), ¿1(1) — невязка на первой итерации (к = 1) для т1 = 1. Ограничение 0,1 ^ т0 ^ 1.

Их применение позволяет сократить число итераций при решении конкретных задач, что видно из сравнения результатов, представленных в следующем разделе в табл. 1 для уравнения Шрёдингера с потенциалом Морзе и табл. 2 при решении уравнения Лежандра.

3. Численные результаты

Программа БЫРМ проверена на ряде тестов [8]. Здесь приведены результаты её применения для решения трёх задач.

Пример 1. Уравнение Шрёдингера с потенциалом Морзе [2]. В этом примере задача (1)-(4) имеет коэффициенты уравнения р(х) = 0, г(х) = 1,

д(х) = -2МВ (е-2«1(«о) - 2е-аг( на интервале (а,Ь) и граничные условия

6,1 = (12 = 1; /1 = ^А - V2МБеа1Х°, /2 = ^А - ^2МБе-а1(ь-х°\ где значения констант а1 = 0.67, х0 = 2.15, М = 4.69, И = 0.1005.

Таблица 1

Хн = 0,435312 Хк/2 = 0,435312 Хк/4 = 0,435312 СТА ~ 4

X Ук Ук/2 Ук/4 ау (25)

0 0,188079е - 1 0,188090е - 1 0,188092е - 1 3, 996

4 0,461810 0,461811 0,461811 3, 999

8 0,492263е - 1 0,492265е - 1 0,492266е - 1 4, 0

12 0, 3613786е - 2 0, 361380е - 2 0,361380е - 2 4, 0

16 0, 258596е - 3 0, 258597е - 3 0,258597е - 3 4, 001

20 0,184721е - 4 0,184723е - 4 0,184723е - 4 4, 004

¿о 60,097359 60,172738 62,171030

го = 0,1 к = 13

¿13 1,23е - 11 1,35е - 11 1, 37е - 11

го = 1 (23) 3

¿3 2,604е - 8 1,551е - 8 1,362е - 8

то = 0,4499 (24) к = 5

¿5 2,877е - 13 2,671е - 13 3,063е - 13

Таблица 2

Хн = -6, 00 Хк/2 = -6, 00 Хк/4 = -6, 00 Ха = -6, 0

X Ук Ук/2 Ук/4 а(25)

-1, 0 1,581136 1., 581138 1,581138 4, 0

-0, 6 0, 632454е - 1 0,632455е - 1 0, 632455е - 1 4, 0

-0,4 -0,411095 -0,411096 -0,411096 4, 001

0,0 -0, 790568 -0, 790569 -0, 790569 4, 0

¿о 39, 847 79,521 158,869

¿13 то = 0,1 8,85е - 11 к = 13 1, 78е - 10 3, 57е - 10

¿4 го = 1 (23) 2,63е - 9 к = 4 5, 26е - 9 1,05е - 8

¿2 то = 0,45 (24) 2,63е - 9 к = 4 9, 01е - 13 2, 07е - 11

Аналитическое решение с безузловой волновой функцией (число узлов п = 0) имеет вид:

Л = 2МБ

где

1_ а1

2л/2МО_

у(х) = С^е-

£ = 25е а1(х Ха), 5 =-, р =-.

а1 а1

Для указанного набора констант а\, хо, М, Б аналитическое значение Ла = 0,4353114734.

С помощью процедуры NEWTON0 получается начальное приближения

дающее невязку 50 < е, е = 0, 0001, и выход из процедуры БЫРМ происходит сразу же после (0) итераций. Для проверки области сходимости метода и достоинств новых алгоритмов выбора то начальное приближение Хо специально сделано плохим: уо = (Сг * у0е^ + С2; Сг = 1; С2 = 0,3); Хо = 1000,435. В

2

первой части табл. 1 приведены вычисленные от этого начального приближения собственные функции уи(х) и собственные значения Ли на последовательности Аи = (к,к/2, к/4) вдвое сгущающихся сеток при значениях а = 0,Ь = 20, к = 0,008333(^ = 2401). Далее показаны начальное значение невязки 50 ~ 60, эмпирически заданное начальное значение параметра то = 0,1, количество итераций к =13 и значения 5 к после окончания итерационного процесса по условию (17). Параметр Тк определяется по формуле (21). В последнем столбце приводятся значения отношения Рунге для г = (\,у(х)):

¿И - ¿И/2 = -'-.

*Н/2 -

Величина выражения (25) должна быть порядка 0(2Р), где р — точность дискретной аппроксимации задачи (1)—(4). Видно, что точность полученного разностного решения задачи (1)—(4) имеет порядок 0(к2).

Во второй части табл. 1 представлены результаты двух расчётов, отличающихся новым способом счета начального значения то по формулам (23) и (24). Видно, что использование алгоритмов (23), (24) привело к сокращению числа итераций с 13 для т0 = 0,1 до 3 и 5, соответственно.

На рис. 1 показан ход итерационного процесса.

0 5 10 15 20

Рис. 1. Аналитика — H4, Аа = 0.4353114, к = 0 — H1, Л0 = 1000.435, к = 8 - H2, As = 270.779, к = 13 — H3, Ахз = 0.4353125, (H3 и H4 — совпадают), к — номер итерации то = 0.1, т^ определяется по формуле (21)

(25)

Пример 2. Уравнение Лежандра. В этом примере р(х) = 0. Задача (1)-(4) имеет вид

2х \

у" (х) - --2 у'(х) - --2У(х) = 0 х G (-1, 1),

1 — X2 1 — X2

граничные условия

У'(х) - f У(х) = 0, х = -1;

у' (ж) + f У(х) = 0, ж = +1.

Аналитические собственные значения данной задачи определяются по формуле Хп = -п(п + 1), а собственные функции могут быть найдены, например, по следующей аналитической формуле в системе MAPLE: рп := LegendreP [п,х\. Ро = 1; р\ = х; р2 = ^(3ж2 - 1),....

В первой части табл. 2 приведены результаты вычисления решения для п = 2, полученные на последовательности А^, = (h,h/2,h/4) вдвое сгущающихся сеток

к = 0, 0125(^ = 1601), подтверждающие порядок 0(Ь?) сходимости приближенного решения задачи (1)—(4). В качестве начального приближения взято специально ухудшенное известное аналитическое решение, которое можно получить также с

помощью процедуры NEWTONO: у0 = (C1 *yJNewt + С2; Ci = 2; С2 = 0,1); Aq = -5.

Показаны начальные значения невязки ¿о.

Во второй части табл. 2 представлены результаты трёх расчётов, отличающихся способом счета начального значения то: то =0,1 и по формулам (23), (24). Здесь к указывает количество итераций, потребовавшихся для вычисления требуемых решений с заданной точностью е = 0,0001. Видно, что использование алгоритмов (23), (24) привело к сокращению числа итераций с 13 для т0 = 0,1 до 4 для обоих новых алгоритмов. Значения параметра и невязки определяются по формулам (21) и (17).

На рис. 2 показан ход итерационного процесса.

Рис. 2. Аналитика — И4, Аа = -6, к = 0 — И1 А0 = -5, к = 4 — И2 А4 = -5, 381, к = 13 — И3, Ахз = -6, 0 (И3 и И4 — совпадают) к — номер итерации то = 0,1 тд. определяется по формуле (21)

Пример 3. Уравнение Шрёдингера

с краевыми условиями у(0) = 0; = 0. Использование процедуры SLIPM от начального приближения Aq = 0, 9 на интервале (0, 20) даёт собственное значение А = 1. Соответствующая волновая функция выдаётся в дискретном виде на сетке по ж. Отметим, что использование процедуры dsolve в СКА Maple для решения этого уравнения даёт общее решение в аналитическом виде, зависящем от А, через функции Уиттекера:

Н1

Н2

НЭ, Н4

(26)

у (ж) = Ci WhittakerM

(27)

Подстановка значения А = 1 в dsolve даёт аналитическое решение у(х) = же х.

4. Задача Ш^турма—Лиувилля в СКА Maple

Для решения задачи Штурма-Лиувилля для ряда уравнений (Лежандра, Уит-текера, Эйлера, Чебышева и др.) [9] в системе Maple при использовании процедуры dsolve можно получить решение в аналитическом виде, зависящем от А. Если мы знаем собственное значение, то с помощью dsolve можно найти соответствующую собственную функцию. Однако, если не известно собственное значение, что обычно бывает на практике при исследовании многих уравнений, то этот подход не работает. Поэтому представляет интерес использование тандема: сначала с помощью SLIPM найти собственное значение, затем использовать его в краевых условиях dsolve для вычисления собственной функции в аналитическом виде, как это показано в третьем примере для уравнения Уиттекера.

5. Заключение

— На базе НАМН создана программа SLIPM на языке MAPLE, предназначенная для численного решения частичной проблемы Штурма-Лиувилля для линейного дифференциального оператора второго порядка с однородными граничными условиями.

— Предложены и программно реализованы два новых алгоритма вычисления начального значения итерационного параметра tq. На приведённых примерах продемонстрированы их достоинства: расширилась область применимости НАМН и сократилось количество итераций.

— На ряде тестов продемонстрирована заданная точность порядка 0(h2) итерационных ньютоновских схем (h — шаг дискретной сетки по Xi). Использование экстраполяции по Ричардсону с использованием результатов на последовательности трёх сгущающихся сеток позволяет повысить точность до порядка 0(h4).

— Возможности графических пакетов СКА MAPLE позволяют организовать хорошую визуализацию итерационного процесса в SLIPM.

— Использование вычисленного с помощью процедуры SLIPM собственного значения в процедуре dsolve открывает новые возможности для получения собственной функции в аналитическом виде.

Литература

1. Puzynin I. V., Puzynina T. P. SLIP1 — Program for the Numerical Solution of the Sturm-Liouville Problem Basing on the Continuous Analog of the Newton Method // Collection of scientific papers in collaboration of JINR, Dubna, USSR and Central Research Institute for Physics, Budapest, Hungary, KFKI-74-34. — 1974. — Pp. 93-112.

2. Puzynin I. V., Puzynina T. P., Strizh T. A. Program for Numerical Solution of the Sturm-Liouville Problem. — JINR Report, P11-87-332. — 1987.

3. Пузынин И. В. и др. Обобщённый непрерывный аналог метода Ньютона для численного исследования некоторых квантово-полевых моделей // ЭЧАЯ. — 1999. — Т. 30. — С. 210-265.

4. Пузынин И. В. и др. О методах вычислительной физики для исследования моделей сложных физических систем // ЭЧАЯ. — 2007. — Т. 38. — С. 144-232.

5. Уилкинсон Дж. Х. Алгебраическая проблема собственных значений. — M.: Наука, 1970.

6. Акишин П. Г., Пузынин И. В. Реализация метода Ньютона в разностной задаче Штурма-Лиувилля. — 1977.

7. Ермаков В. В., Калиткин Н. Н. Оптимальный шаг и регуляризация метода Ньютона // ЖВМиМФ. — 1981. — Т. 21. — С. 491.

8. Bailey P. B, Everitt W. N., Zettl A. The SLEIGN2 Sturm-Liouville Code // ACM Trans. Math. Software. — 2001. — Vol. 21. — Pp. 143-192.

98

ny3MHHH H. B., ny3MHHHa T. n., TxaK B. H.

9. Lopez R. J. Emeritus Professor of Mathematics and Maple Fellow Classroom Tips and Techniques Eigenvalue Problems for ODEs. Part1-Part3. — Maplesoft.

UDC 517.9, 519.6, 539.12

SLIPM — MAPLE Programm for Numerical Solution of Sturm—Liouville Partial Problem with the Help of the Continuous Analogue of Newton's Method

I. V. Puzynin, T. P. Puzynina, V. T. Thach

A computer program SLIPM (Sturm-Liouville Problem in MAPLE) has been designed in the language of computer algebras MAPLE. It is intended for the numerical solution of Sturm-Liouville partial problem with the help of the continuous analogue of Newton's method, i.e., for calculating the eigenvalue of a linear second-order differential operator and a corresponding eigenfunction satisfying homogeneous boundary conditions. SLIPM is the development of the computer codes SLIP1 and SLIPH4 written in the Fortran language; it is added by a new procedure of calculating an initial approach to the solution and by new ways of calculating the initial value of the iterative parameter t0.

Key words and phrases: computational mathematics, continuous analogue of Newton's method, computer algebras system, MAPLE.