Numerical methods for solving direct and inverse Sturm–Liouville problems

Full Text

Постановка задач

Будет рассматриваться задача Штурма–Лиувилля второго порядка с краевыми условиями Дирихле, приведённая в [1]:

( -y^{\prime\prime}+q\left(x\right)y=\lambda y,x\in\left(0,\pi\right),\left(1.1\right) )

( y(0)=0,y(π)=0,(1.2) )

где ( \lambda ) – спектральный параметр, ( q\left(x\right)\in L_2\left(0,\pi\right) ) – вещественный потенциал.

Интервал ( \left(0,\pi\right) ) и краевые условия Дирихле были взяты для определённости, так как к данному интервалу можно привести произвольный интервал путём линейной замены, а краевые условия принципиально не влияют на принцип решения задач.

Собственные значения краевой задачи – это такие значения спектрального параметра ( \lambda ), при которых решения ( y ) будут нетривиальными. Такие решения называются собственными функциями. У краевой задачи Штурма–Лиувилля существует множество собственных чисел ( \{{{\lambda}_n}\}_{n=1}^\infty ), называемое спектром. В работе [1] для каждого собственного значения определены весовые числа

( \alpha_n=\int_{0}^{\pi}y_n^2dx,\left(1.3\right) )

где ( y_n ) – собственная функция при ( \lambda_n ).

Множество ( \{{\lambda}_n,{\alpha_n}\}_{n=1}^\infty ) называется спектральными данными задачи (1.1)–(1.2). Постановим две задачи для (1.1)–(1.2).

Прямая задача: для (1.1)–(1.2) дан потенциал ( q(x) ). Необходимо найти ( \{{\lambda}_n,{\alpha_n}\}_{n=1}^N ).

Обратная задача: положим ( B(q) ) – задача (1.1)–(1.2) с известными спектральными данными ( \{{\lambda}_n,{\alpha_n}\}_{n=1}^N ) и неизвестным потенциалом ( q(x) ). Необходимо построить ( q(x) ).

Прямая задача возникает при решении уравнений математической физики методом разделения переменных, описанным в [2]. А обратная задача позволяет определить неизвестные свойства среды по наблюдаемым величинам. В работе [3] приведены примеры приложений. Например, с помощью обратной задачи можно распознать тип крепления мембраны по собственным частотам радиальных колебаний. Другой пример – акустическая диагностика вала с двумя концевыми дисками. Она позволяет выявить не только неисправности, но и обнаружить только развивающиеся дефекты, и определить необходимость ремонта, его объём и сроки.

Точное аналитическое решение прямой и обратной задачи можно получить только в частных случаях. В общих случаях используются численные методы или асимптотические формулы для приближенного решения. Некоторые численные методы решения задач предложены в работах [4;5;6]. В статье [4] рассматривается метод групп Ли для решения прямых задач Штурма-Лиувилля высших порядков. В [6] разработан численный алгоритм решения обратной задачи на основе метода спектральных отображений, который является наиболее эффективным из теоретических подходов в теории обратных спектральных задач (см. [7]). Однако в [6] рассмотрена задача с краевыми условиями Неймана. В данной статье предложенный подход реализован для случая условий Дирихле, который обладает рядом технических отличий, а также более детально проведен численный анализ. Для прямой и обратной задач были составлены и реализованы алгоритмы решения, проведены эксперименты на тестовых примерах, доказывающие корректность работы данных методов.

Решение прямой задачи

Для составления алгоритма использовались статья [1]. Определим ( S(x, \lambda) ) как решение задачи (1.1) с начальными условиями

 ( S(0, \lambda) = 0, S\prime(0, \lambda)=1. (2.1) )

Первое начальное условие соответствует первому краевому условию (1.2). Второе начальное условие взято для определённости.

Определим характеристическую функцию

 ( \Delta\left(\lambda\right)=S\left(\pi,\lambda\right). )

Если ( \Delta\left(\lambda\right)=0 ), то решение ( S\left(x,\lambda\right) ) удовлетворяет краевым условиям (1.2), и, следовательно, ( \lambda ) и ( S(x, \lambda) ) – собственное значение и собственная функция задачи (1.1)–(1.2) соответственно. Для вычисления решения задачи Коши (1.1)–(2.1) использовался классический метод Рунге–Кутта четвёртого порядка, приведённый в [8]. Для использования этого метода уравнение второго порядка (1.1) было сведено к системе уравнений первого порядка

( \left\{ \begin{aligned}
u_1^\prime=u_2 \\
u_2^\prime=q\left(x\right)u_1-\lambda u_1
\end{aligned} \right. )

где ( u_1 =y, u_2 = y\prime ).

Начальные условия (2.1) примут вид

( u_1\left(0\right)=0,\ \ u_2\left(0\right)=1. )

Таким образом, задача свелась к поиску корней уравнения

( \Delta\left(\lambda\right)=0. )

Их поиск осуществлялся методом дихотомии.

Пусть задана точность для метода дихотомии ( \epsilon ). Для обеспечения заданной точности в методе дихотомии необходимо подобрать точность ( \epsilon^\prime ) для метода Рунге–Кутта. Разложим ( S\left(\pi,\lambda\right) ) в ряд Тейлора в окрестности ( \lambda_n )

( S\left(\pi,\lambda\right)=S\left(\pi,\lambda_n\right)+\dot{S}\left(\pi,\lambda_n\right)\left(\lambda-\lambda_n\right)+O\left(\left(\lambda-\lambda_n\right)^2\right). )

Так как ( S\left(\pi,\lambda_n\right)=0 ):

( S\left(\pi,\lambda\right)\approx\dot{S}\left(\pi,\lambda_n\right)\left(\lambda-\lambda_n\right), )

( \lambda-\lambda_n\approx\frac{S\left(\pi,\lambda\right)}{\dot{S}\left(\pi,\lambda_n\right)}. )

В силу того, что ( \left|S\left(\pi,\lambda\right)-S\left(\pi,\lambda_n\right)\right|<\epsilon^\prime ), ( S\left(\pi,\lambda_n\right)=0 ), и ( \left|\lambda-\lambda_n\right|<\epsilon ),  получим

( \epsilon\sim\frac{\epsilon^\prime}{\frac{d}{d\lambda}S\left(\pi,\lambda_n\right)}\ \ \ \rightarrow\ \ \ \epsilon^\prime\sim\ \epsilon\frac{d}{d\lambda}S\left(\pi,\lambda_n\right). )

( \epsilon^\prime\sim\epsilon\frac{d}{d\lambda}\Delta\left(\lambda_n\right). )

Так как значения ( S\left(\pi,\lambda_n\right) ) и ( \lambda_n ) неизвестны, будут использоваться асимптотические формулы, приведённые в [7]:

( S\left(\pi,\lambda\right)\approx\frac{sin\left(\rho\pi\right)}{\rho}, )

( ρ_n≈n,(2.2) )

где ( \rho=\sqrt\lambda ).

Тогда положим

( \epsilon^\prime=\epsilon\frac{d}{d\lambda}\left(\frac{sin\left(\sqrt{\lambda_n}\pi\right)}{\sqrt{\lambda_n}}\right)=\epsilon\frac{\pi\sqrt{\lambda_n}cos\left(\pi\sqrt{\lambda_n}\right)-sin\left(\pi\sqrt{\lambda_n}\right)}{2\lambda_n\sqrt{\lambda_n}}=\epsilon\frac{\pi}{2n^2}. )

Так как ( \frac{\pi}{2n^2}\ <\ 1 ), то ( \epsilon^\prime<\epsilon ).

Экспериментальным путем было установлено, что ( h_x\rightarrow0 ) при ( \lambda\rightarrow\infty ) и при фиксированном ( \epsilon ). Для ( \{{\lambda_n}\}_{n=1}^{20} ) шаг ( h_x ) для метода Рунге–Кутта вычисляется методом Рунге при ( \lambda_{20} ). Так как оно неизвестно, в силу (2.2) положим ( \lambda_{20}={20}^2 ). Далее, при ( n\vdots\ 20 ) шаг ( h_x ) уменьшается в 2 раза. На основе данных рассуждений был составлен алгоритм решения прямой задачи.

Для вычисления интеграла в формуле (1.3) использовалась формула Ньютона из [8]:

( \int_{a}^{b}f\left(x\right)dx\approx\frac{3h}{8}\left(y_0+y_{3m}+2\left(y_3+\ldots+y_{3m-3}\right)+3\left(y_1+y_2+\ldots y_{3m-2}+y_{3m-1}\right)\right),\left(2.3\right) )

где ( h=\frac{b-a}{n}=\frac{b-a}{3m}. )

Собственные значения могут быть отрицательными. Чтобы их учесть, можно воспользоваться следующими фактами:

1. ( \lambda_n\geq0 ) при ( q\left(x\right)\geq0 ).
2. При прибавлении константы ( c ) к потенциалу ( q ), собственные числа изменятся на ( c ):

( -y^{\prime\prime}+\left(q\left(x\right)+c\right)y=\left(\lambda+c\right)y. )

Если потенциал ( q(x) ) допускает отрицательные значения, то найдем минимальное значение потенциала ( q_{min} ) и составим задачу с неотрицательным потенциалом

( -y^{\prime\prime}+\left(q\left(x\right)-q_{min}\right)y=\left(\lambda-q_{min}\right)y.\left(2.4\right) )

Получив собственные значения ( {\overline{\lambda}}_n ) для задачи (2.4), можно получить собственные значения для исходной задачи по формуле ( \lambda_n={\overline{\lambda}}_n+q_{min} ).

Алгоритм 1. Положим, что задан потенциал ( q(x) ), N для множества ( \{{\lambda}_n,{\alpha_n}\}_{n=1}^N ) и точность ( \epsilon ).

1. Вычислить ( \epsilon^\prime ). Определить ( i=0 ), вычислить шаг ( h_x ) методом Рунге при ( \lambda_{20} ) и вычислить ( \Delta\left(\lambda_0\right) ). Определить ( q_{min} ). Перейти к шагу 2.
2. Если ( q_{min}<0 ), то положить ( c=q_{min} ), иначе ( c=0 ). Определить потенциал ( q\left(x\right)=q\left(x\right)-c. ). Перейти к шагу 3.
3. Положить ( i=i+1 ) и вычислить ( \Delta(\lambda_i) ). Если ( sign\left(\Delta\left(\lambda_{i-1}\right)\right)\neq sign\left(\Delta\left(\lambda_i\right)\right) ), то перейти к шагу 4. Иначе выполнить шаг 3.
4. Методом дихотомии уточнить корень уравнения на интервале ( \left[\lambda_{i-1},\lambda_i\right] ) до необходимой точности. Положить .( \lambda_n=\lambda_n+c ) Перейти к шагу 5.
5. Вычислить весовое число ( \alpha_i ) используя формулу (2.3). Если ( n\vdots\ 20, ) то ( h_x=\frac{h_x}{2} ). Перейти к шагу 3, пока не будет найдено ( \lambda_N ).

Численные эксперименты

Данный алгоритм был реализован на С++. Для проверки корректности алгоритма проводились численные эксперименты при ( N=10,\ \ \epsilon={10}^{-4} ). Для первого эксперимента был взят частный случай при ( q\left(x\right)=0 ). В данном случае прямую задачу можно решить аналитически и получить формулы для сравнения с численными результатами (см. табл. 1):

( \lambda_n=n^2,\ \ n\in\mathbb{N}, )

( \alpha_n=\int_{0}^{\pi}{\frac{\sin^2{\left(\sqrt{\lambda_n}x\right)\ }}{\lambda_n}dx}=\frac{\pi}{2n^2}. )

Таблица 1

Решение прямой задачи при ( q(x) = 0 )

n	( {\overline{\lambda_n}} )	( \lambda_n )	( \left|\overline{\lambda_n} - \lambda_n \right| )	( \overline{\alpha_n} )	( \alpha_n )	( \left| \overline{\alpha_n} - \alpha_n \right| )
1	1,0000305	1	3,0518( \cdot 10^{-5} )	1,5707229	1,5707963	7,34( \cdot 10^{-5} )
2	4,0000305	4	3,0518( \cdot 10^{-5} )	0,3926931	0,3926991	0,6( \cdot 10^{-5} )
3	9,0000305	9	3,0518( \cdot 10^{-5} )	0,1745307	0,1745329	0,022( \cdot 10^{-5} )
4	16,0000305	16	3,0518( \cdot 10^{-5} )	0,0981730	0,0981748	0,018( \cdot 10^{-5} )
5	25,0000305	25	3,0518( \cdot 10^{-5} )	0,0628298	0,0628319	0,021( \cdot 10^{-5} )
6	36,0000305	36	3,0518( \cdot 10^{-5} )	0,0436318	0,0436332	0,014( \cdot 10^{-5} )
7	49,0000305	49	3,0518( \cdot 10^{-5} )	0,0320551	0,0320571	0,020( \cdot 10^{-5} )
8	64,0000305	64	3,0518( \cdot 10^{-5} )	0,0245422	0,0245437	0,015( \cdot 10^{-5} )
9	81,0000305	81	3,0518( \cdot 10^{-5} )	0,0193902	0,0193925	0,023( \cdot 10^{-5} )
10	100,0000305	100	3,0518( \cdot 10^{-5} )	0,0157060	0,0157080	0,020( \cdot 10^{-5} )

Здесь ( {\overline{\lambda}}_n ) и ( {\overline{\alpha}}_n ) – собственные значения и весовые числа, полученные численным методом, ( \lambda_n ) и ( \alpha_n ) – собственные значения и весовые числа, полученные аналитически.

Второй эксперимент проводился при ( q\left(x\right)=-2sin\left(3x\right) ). Численные результаты сравнивались с асимптотической формулой (табл. 2) из работы [7]:

( \sqrt{\lambda_n}=n+\frac{\omega}{\pi n}+o\left(\frac{1}{n}\right),\left(2.5\right) )

где ( \omega=\frac{1}{2}\int_{0}^{\pi}q\left(s\right)ds ).

При ( q\left(x\right)=-2sin\left(3x\right) ) получим асимптотику

( \lambda_n\approx\left(n-\frac{2}{3\pi n}\right)^2. )

Таблица 2

 Решение прямой задачи при ( q(x) = -2sin(3x) )

n	( \overline{\lambda}_n )	( \lambda_n )	( \left| \overline{\lambda}_n - \lambda_n \right| )	( {\overline{{\alpha}}}_{n} )
1	1,102203	0,620618	0,481585	0,7778649
2	2,958832	3,586844	0,628012	0,3323804
3	8,613129	8,580590	0,032539	0,1808469
4	15,531403	15,578401	0,046998	0,0988727
5	24,553864	24,577388	0,023524	0,0633085
6	35,560517	35,576837	0,01632	0,0438726
7	48,564728	48,576505	0,011777	0,0321888
8	63,567413	63,576290	0,008877	0,0246216
9	80,569183	80,576142	0,006959	0,0194401
10	99,570465	99,576037	0,005572	0,0157389

На основе этих двух численных экспериментов можно сделать вывод, что алгоритм работает корректно.

Обратная задача

В данном разделе использовался подход, предложенный в работе [6], который основан на методе спектральных отображений. Предложенный алгоритм решения обратной задачи был изменен для разработки программы.

Пусть дана задача ( B=B(q) ). Для построения ( q(x) ) необходимо выбрать задачу ( \widetilde{B}=B\left(\widetilde{q}\right) ). Для нее необходимо подобрать такой ( \widetilde{q} ), чтобы ( \omega = \widetilde{\omega} ). Так как потенциал ( q(x) ) неизвестен для вычисления ( \omega ), можно получить приближённое значение из асимптотики (2.5)

( \omega\approx\pi n\left(\sqrt{\lambda_n}-n\right). )

Идея метода заключается в итерационной аппроксимации потенциала в точке ( x ) по формуле

( q_r(x)=\widetilde{q}(x)-2ϵ_N'(x),\left(3.1\right) )

где ( r ) – текущая итерация аппроксимации,

( \epsilon_N\left(x\right)=\sum_{n=0}^{N}\left(\frac{1}{\alpha_n}S\left(x_i,\lambda_n\right)\widetilde{S}\left(x_i,\lambda_n\right)-\frac{1}{{\widetilde{\alpha}}_n}S\left(x_i,{\widetilde{\lambda}}_n\right)\widetilde{S}\left(x_i,{\widetilde{\lambda}}_n\right)\right). )

Следовательно,

( \epsilon_N^\prime\left(x\right)=\sum_{n=0}^{N}\left({\frac{1}{\alpha_n}\left(S\left(x_i,\lambda_n\right)\widetilde{S}\left(x_i,\lambda_n\right)\right)}^\prime-{\frac{1}{{\widetilde{\alpha}}_n}\left(S\left(x_i,{\widetilde{\lambda}}_n\right)\widetilde{S}\left(x_i,{\widetilde{\lambda}}_n\right)\right)}^\prime\right). )

Далее, с помощью свойства вычисления производной произведения можно получить производные слагаемых членов суммы.

Отметим, что для вычисления ( \epsilon_N^\prime ) производные ( S^\prime\left(x,\lambda\right) ) вычисляются в методе Рунге–Кутта и нет необходимости вычислять производные отдельно.

Используя ( \widetilde{B} ), можно использовать следующий алгоритм для решения обратной задачи.

Алгоритм 2. Положим, даны спектральные данные ( \{{\lambda}_n,{\alpha_n}\}_{n=1}^N ), точность ( \epsilon ) и количество узлов ( m ).

1. Подобрать ( \widetilde{B} ) и вычислить ( \{\widetilde{\lambda_n},{\widetilde{\alpha_n}\}}_{n=1}^N ). Перейти к шагу 1.
2. Положить начальное приближение ( q_0\left(x_i\right)=\widetilde{q}\left(x_i\right) ) и номер итерации приближения ( r=1. ) Перейти к шагу 3.
3. Вычислить методом Рунге–Кутта решения ( S\left(x_i,\lambda_n\right), \widetilde{S}\left(x_i,\ \lambda_n\right), S\left(x_i,{\widetilde{\lambda}}_n\right), \widetilde{S}\left(x_i,{\widetilde{\lambda}}_n\right),i=\overline{0,m},n=\overline{1,N} ) и их производные. Перейти к шагу 4.
4. Вычислить ( \epsilon_N^\prime\left(x_i\right),\ i=\overline{0,\ m} ) по формуле (3.2). Перейти к шагу 5.
5. Получить аппроксимацию на текущей итерации по формуле (3.1). Если ( max\left|q_r\left(x_i\right)-q_{r-1}\left(x_i\right)\right|<\epsilon,i=\overline{0,m} ), то закончить вычисления. Иначе положить ( r=r+1 ) и перейти к шагу 3.

Численные эксперименты

Данный алгоритм также был реализован на С++. Эксперимент заключался в задании ( q(x) ) и ( \widetilde{q}(x) ), решении для них прямой задачи, и решении обратной задачи по полученным спектральным данным. Все эксперименты проводились при ( N=10, m=10, \epsilon={10}^{-4}. )

Первый эксперимент проводился для ( q\left(x\right)=cos\left(x\right) ). Было подобрано начальное приближение ( \widetilde{q}\left(x\right)=1-\frac{2x}{\pi} ). Полученные результаты сравнивались с исходным заданным потенциалом (см. табл. 3). Во время вычислений было произведено 6 итераций.

Таблица 3

Решение обратной задачи при ( q(x) = cos(x) )

( x )	( cos(x) )	( q(x) )	( \left| cos(x) - q(x) \right| )
0,0000	1,0000	1,0000	0,0000
0,3491	0,93969	0,93928	4,1756( \ \cdot{10}^{-4} )
0,6981	0,76604	0,76587	1,6992( \ \cdot{10}^{-4} )
1,0472	0,50000	0,50002	2,3917( \ \cdot{10}^{-5} )
1,3963	0,17365	0,17407	4,2143( \ \cdot{10}^{-4} )
1,7453	-0,17365	-0,17320	4,4679( \ \cdot{10}^{-4} )
2,0944	-0,50000	-0,50014	1,3881( \ \cdot{10}^{-4} )
2,4435	-0,76604	-0,76686	8,1907( \ \cdot{10}^{-4} )
2,7925	-0,93969	-0,94107	1,3792( \ \cdot{10}^{-3} )
3,1416	-1,0000	-0,98439	1,5613( \ \cdot{10}^{-2} )

Аналогично была решена обратная задача при ( q(x)=2x ) с начальным приближением ( \widetilde{q}\left(x\right)=\frac{x^2}{\pi}+\frac{4x}{3} ) (см. табл. 4). Для достижения необходимой точности потребовалось 8 итераций.

Таблица 4

Решение обратной задачи при ( q(x) = 2x )

( x )	( 2x )	( q(x) )	( |2x - q(x)| )
0,0000	0,0000	0,00000	0,0000
0,3491	0,698132	6,98809	6,7692( \cdot{10}^{-4} )
0,6981	1,396263	1,39696	7,0146( \cdot{10}^{-4} )
1,0472	2,094395	2,09428	1,1679( \cdot{10}^{-4} )
1,3963	2,792527	2,79133	1,1951( \cdot{10}^{-3} )
1,7453	3,490659	3,48943	1,2264( \cdot{10}^{-3} )
2,0944	4,188790	4,18948	6,9433( \cdot{10}^{-4} )
2,4435	4,886922	4,89006	3,1411( \cdot{10}^{-3} )
2,7925	5,585054	5,58957	4,5152( \cdot{10}^{-3} )
3,1416	6,283185	6,22588	5,7307( \cdot{10}^{-2} )

Заметим, что на правом конце интервала точность снижается. Для повышения точности необходимо подбирать потенциалы, удовлетворяющие краевым условиям ( q\left(0\right)=0, q\left(\pi\right)=0 ).

Заключение

В данной работе был сформулирован алгоритм для решения прямой задачи на основе работ [1; 7], разработана программа и проведены численные эксперименты для проверки корректности алгоритма. Также был изучен метод спектральных отображении и метод решения обратной задачи, предложенный в [6], разработана программа и проведены численные эксперименты решения обратной задачи. В ходе разработки программ были реализованы методы Рунге–Кутта и вычисления интеграла по формуле Ньютона, метод Рунге для вычисления шага по  для методов Рунге–Кутта и вычисления интеграла. Разработанные программы можно использовать для решения прямой и обратной задач Штурма–Лиувилля второго порядка. В дальнейшем на основе полученных методов и алгоритмов могут быть разработаны подходы решения прямых и обратных задач Штурма–Лиувилля высших порядков.

About the authors

Rafael Abeldinov

Samara National Research University named after Academician S.P. Korolev

Author for correspondence.
Email: raf2003@bk.ru

Студент

Russian Federation, 443086, Russia, Samara, Moskovskoye shosse, 34.

References

Supplementary files