Решение краевых ОДУ

 х0 — левая граница расчетного интервала.

x1 — правая граница  расчетного интервала. 

 xf — точка внутри интервала.

Dо(х,у) — векторная функция, описывающая систему N ОДУ размера Nx1 и двух аргументов — скалярного х и векторного у. При этом у — это неизвестная векторная функция аргумента х того же размера Nx1.

load1(x0,z) — векторная  функция размера Nx1 левых граничных  условий, причем недостающие начальные  условия поименовываются соответствующими компонентами векторного аргумента z.

load2(x1,z) — векторная  функция размера Nx1 левых граничных  условий, причем недостающие начальные  условия поименовываются соответствующими компонентами векторного аргумента z.

 score (xf, у) — векторная функция размера Nx1, выражающая внутреннее условие для векторной функции у в точке xf.  
 

Обычно функция  bvaifit применяется для задач, в которых производная решения имеет разрыв во внутренней точке xf. Некоторые из таких задач не удается решить обычным методом пристрелки, поэтому приходится вести пристрелку сразу из обеих граничных точек. В этом случае дополнительное внутреннее условие в точке xf является просто условием сшивки в ней левого и правого решений. Поэтому для данной постановки задачи оно записывается В форме score (xf, у) :=у.  

Рассмотрим действие функции bvaifit на знакомом примере модели взаимодействия пучков света (см. рис. 10.1), предположив, что в промежутке между xf=0.5 и x1=1.0 находится другая, оптически более плотная среда с другим коэффициентом ослабления излучения а(х)=3. Соответствующая краевая задача решена в листинге 10.3, причем разрывный показатель ослабления определяется в его второй строке.  

Листинг 10.3. Краевая  задача с граничным условием в  промежуточной точке

 

 

Система уравнений  и левое краевое условие вводится так же, как и в предыдущем листинге для функции sbval. Обратите внимание, что таким же образом записано и правое краевое условие. Для того чтобы ввести условие отражения на правой границе, пришлось определить еще один неизвестный пристрелочный параметр z20. Строка листинга, в которой определена функция score, задает условие стрельбы — сшивку двух решений в точке xf. В самой последней строке листинга выдан ответ — определенные численным методом значения обоих пристрелочных параметров, которые объединены в вектор а (мы применили в предпоследней строке операцию транспонирования, чтобы результат получился в форме вектора, а не матрицы-строки). Для корректного построения графика решения лучше составить его из двух частей — решения задачи Коши на интервале (x0,xf) и другой задачи Коши для интервала (xf,x1). Реализация этого способа приведена в листинге 10.4, который является продолжением листинга 10.3. В последней строке листинга 10.4 выведено значение второй искомой функции на правой границе интервала. Всегда полезно проконтролировать, что оно совпадает с соответствующим пристрелочным параметром (выведенным в последней строке листинга 10.3).  

Листинг 10.4. Решение  краевой задачи (продолжение листинга 10.3)

   

Решение краевой  задачи приведено на рис. 10.5. С физической точки зрения естественно, что интенсивность  света уменьшается быстрее по мере распространения в более  плотной среде в правой половине расчетного интервала. В средней  точке xf=0.5, как и ожидалось, производные  обоих решений имеют разрыв.  
 

ПРИМЕЧАНИЕ  

Еще один пример решения  краевых задач с разрывными коэффициентами ОДУ приведен в справочной системе  Mathcad.  
 
 
 
 

Рис. 10.5. Решение краевой  задачи с разрывом в xf=0.5 (продолжение  листингов 10.3—10.4)

 

Ради справедливости необходимо заметить, что разобранную  краевую задачу легко решить и  с помощью функции sbval, заменив в листинге 10.2 зависимость а (х) на третью строку листинга 10.3. В этом случае листинг 10.2 даст в точности тот же ответ, что показан на рис. 10.5. Однако в определенных случаях (в том числе из соображений быстроты расчетов) удобнее использовать функцию bvaifit, т. е. вести пристрелку с обеих границ интервала.  
 

СОВЕТ  

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

Задачи на собственные  значения — это краевые задачи для системы ОДУ, в которой  правые части зависят от одного или  нескольких параметров X. Значения этих параметров неизвестны, а решение  краевой задачи существует только при определенных λk, которые называются собственными значениями (eigenvalues) задачи. Решения, соответствующие этим λk, называют собственными функциями (eigenflmctions) задачи. Правильная постановка таких задач требует формулировки количества граничных условий, равного сумме числа уравнений и числа собственных значений. Физическими примерами задач на собственные значения являются, например, уравнение колебаний струны, уравнение Шредингера в квантовой механике, уравнения волн в резонаторах и многие другие.  

С вычислительной точки  зрения, задачи на собственные значения очень похожи на рассмотренные выше краевые задачи. В частности, для  многих из них также применим метод  стрельбы (см. разд. 10.1.2). Отличие заключается в пристрелке не только по недостающим левым граничным условиям, но еще и по искомым собственным значениям. В Mathcad для решения задач на собственные значения используются те же функции sbval и bvaifit. В их первый аргумент, т. е. вектор, присваивающий начальные значения недостающим начальным условиям, следует включить и начальное приближение для собственного значения.  

Рассмотрим методику решения на конкретном примере определения  собственных упругих колебаний  струны. Профиль колебаний струны у(х) описывается линейным дифференциальным уравнением второго порядка

 

  

где р(х) и q(x) — жесткость и плотность, которые, вообще говоря, могут меняться вдоль струны. Если струна закреплена на обоих концах, то граничные условия задаются в виде у(0)=у(1)=0. Сформулированная задача является частным случаем задачи Штурма—Лиувилля. Поскольку решается система двух ОДУ, содержащая одно собственное значение x, то по идее задача требует задания трех (2+1) условий. Однако, как легко убедиться, уравнение колебаний струны — линейное и однородное, поэтому в любом случае решение у(х) будет определено с точностью до множителя. Это означает, что производную решения можно задать произвольно, например, у' (0)=1, что и будет третьим условием. Тогда краевую задачу можно решать как задачу Коши, а пристрелку вести только по одному параметру — собственному значению.  

Процедура поиска первого  собственного значения представлена в  листинге 10.5.

Листинг 10.5. Решение  задачи о собственных колебаниях струны

  

В первых двух строках  листинга определяются функции, входящие в задачу, в том числе р' (х) :=0, и границы расчетного интервала (0,1). В третьей строке дается начальное приближение к собственному значению λ0, в четвертой вводится система ОДУ. Обратите внимание, что она состоит не из двух, а из трех уравнений. Первые два из них определяют эквивалентную (10.3) систему ОДУ первого порядка, а третье необходимо для задания собственного значения в виде еще одного компонента у2 искомого вектора у. Поскольку по определению собственное значение постоянно при всех х, то его производная должна быть приравнена нулю, что отражено в последнем уравнении. Важно также, что во втором из уравнений собственное значение записано как у2, поскольку является одним из неизвестных.  

В следующих двух строках листинга задается левое  граничное условие, включающее и  недостающее условие на собственное  значение для третьего уравнения, и  правое граничное условие у0=0. В предпоследней строке листинга обычным образом применяется функция sbval, а в последней выводится результат ее работы вместе с известным аналитически собственным значением п2-я2. Как легко убедиться, мы нашли первое собственное значение для n=1, а чтобы найти другие собственные значения, необходимо задать другие начальные приближения к ним (в третьей строке листинга 10.5). Например, выбор λ0=50 приводит ко второму собственному значению 22π2, а λ0=100 — к третьему 32π2.  

Чтобы построить  график соответствующей собственной  функции, надо добавить в листинг  строку, программирующую решение  задачи Коши, например, такую: u:=rkfixed(load(a,A) ,а,b, 100,D) . Полученные кривые показаны на рис. 10.6 в виде коллажа трех графиков, рассчитанных для трех собственных значений.  
 

ПРИМЕЧАНИЕ  

Примеры решения  нескольких задач на собственные  значения можно найти в различных  Ресурсах Mathcad.  
 
 
 
 

Рис. 10.6. Первые три  собственные функции задачи колебаний  струны (коллаж трех графиков, рассчитанных листингом 10.5 для разных λ0)

 

Многие краевые  задачи не поддаются решению методом  стрельбы. Однако в Mathcad 12 других встроенных алгоритмов нет. Тем не менее это не означает, что по-другому решать краевые задачи невозможно, ведь другие численные алгоритмы несложно запрограммировать самому пользователю. Рассмотрим возможную реализацию наглядного метода, называемого разностным, которым можно решать краевые задачи как для ОДУ, так и для дифференциальных уравнений в частных производных.

Разберем идею разностного  метода решения краевых задач  на примере взаимодействия световых пучков (см. рис. 10.1), переобозначив в системе (10.1) интенсивность излучения вправо на Y, а интенсивность излучения влево на у (просто в целях удобства, чтобы не писать в дальнейшем индекс). Суть метода заключается в покрытии расчетного интервала сеткой из N точек. Тем самым определяются (м-l) шагов (рис. 10.7). Затем надо заменить дифференциальные уравнения исходной краевой задачи аппроксимирующими их уравнениями в конечных разностях, выписав соответствующие разностные уравнения для каждого 1-го шага. В нашем случае достаточно просто заменить первые производные из (10.1) их разностными аналогами (такой метод называется еще методом Эйлера):

 

 

  

Рис. 10.7. Сетка, покрывающая  расчетный интервал  
 
 

ПРИМЕЧАНИЕ  

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

Получилась система (по числу шагов) 2(N-1) разностных линейных алгебраических уравнений с 2N неизвестными Yi и уi. Для того чтобы она имела единственное решение, надо дополнить число уравнений до 2N.  

Это можно сделать, записав в разностном виде оба  граничных условия: \  

Y0=10, YN=RYN. (10.5)  

Сформированная полная система алгебраических уравнений  называется разностной схемой, аппроксимирующей исходную краевую задачу. Обратите внимание, что правые части разностных уравнений системы (10.4) на каждом шаге записаны для левой границы шага. Такие разностные схемы называют явными, т. к. все значения Yi+1 и yi+1 находятся  в левой части уравнений. Полученную явную разностную схему легко  записать в матричной форме  

Аz=В, (10.6)  

где z — неизвестный вектор, получающийся объединением векторов Y и у. Решив систему (10.6), мы получим решение краевой задачи.  
 

ПРИМЕЧАНИЕ  

На самом деле все несколько сложнее, поскольку, вообще говоря, необходимо еще доказать, что, во-первых, разностная схема действительно  аппроксимирует дифференциальные уравнения  и, во-вторых, при N->∞ разностное решение  действительно сходится к дифференциальному.  
 
 

Процесс решения  системы разностных уравнений называют также реализацией разностной схемы. Программа, которая решает рассматриваемую  краевую задачу разностным методом, приведена в листинге 10.6.  

Листинг 10.6. Реализация явной разностной схемы

  

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

В первой строке листинга определяются функции и константы, входящие в модель, во второй задается число точек сетки N=5 и ее равномерный  шаг. Следующие две строки определяют матричные коэффициенты, аппроксимирующие уравнения для Yi; а пятая и шестая — для уi. Седьмая и восьмая строки листинга задают соответственно левое и правое граничное условие, а строки с девятой по одиннадцатую — правые части системы (10.6). В следующей строке завершается построение матрицы А вырезанием из нее левого нулевого столбца. В предпоследней строке листинга применена встроенная функция isolve для решения системы (10.6), а в последней выведены рассчитанные ею неизвестные граничные значения. Графики решения приведены на рис. 10.8, причем первые N элементов итогового вектора есть вычисленное излучение вперед, а последние N элементов — излучение назад.