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

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

Краевые задачи во множестве  практических приложений часто зависят  от некоторого числового параметра. При этом решение существует не для  всех его значений, а лишь для  счетного их числа. Такие задачи называют задачами на собственные значения (см. разд. 10.3).  

Несмотря на то, что, в отличие от задач Коши для  ОДУ, в Mathcad не предусмотрены встроенные функции для решения жестких краевых задач, с ними все-таки можно справиться, применив программирование разностных схем, подходящих для решения задач этого класса (см. разд. 104). О подходе к решению нелинейных краевых задач написано в конце главы (см. разд. 10.5).

Постановка краевых  задач для ОДУ отличается от задач Коши, рассмотренных в главе 9, тем, что граничные условия для них ставятся не в одной начальной точке, а на обеих границах расчетного интервала. Если имеется система N обыкновенных дифференциальных уравнений первого порядка, то часть из N условий может быть поставлена на одной границе интервала, а оставшиеся условия — на противоположной границе.  
 

ПРИМЕЧАНИЕ  

Дифференциальные  уравнения высших порядков можно  свести к эквивалентной системе  ОДУ первого порядка (см. главу 9).  
 
 

Чтобы лучше понять, что из себя представляют краевые задачи, рассмотрим их постановочную часть на конкретном физическом примере модели взаимодействия встречных световых пучков. Предположим, что надо определить распределение интенсивности оптического излучения в пространстве между источником (лазером) и зеркалом, заполненном некоторой средой (рис. 10.1). Будем считать, что от зеркала отражается R-Я часть падающего излучения (т. е. его коэффициент отражения равен R), а среда как поглощает излучение с коэффициентом ослабления а(х), так и рассеивает его. Причем коэффициент рассеяния назад равен г(х). В этом случае закон изменения интенсивности у0(х) излучения, распространяющегося вправо, и интенсивности y1 (х) излучения влево определяется системой двух ОДУ первого порядка:

  

Для правильной постановки задачи требуется, помимо уравнений, задать такое же количество граничных условий. Одно из них будет выражать известную  интенсивность излучения 10, падающего  с левой границы х=0, а второе — закон отражения на его правой границе x=1:  

y0(0) = 10; (10.2)  

y1(l) = Rxy0(l).  
 

Рис. 10.1. Модель краевой  задачи  
 

Полученную задачу называют краевой (boundary value problem), поскольку условия поставлены не на одной, а на обеих границах интервала (0,1).  И, в связи с этим, их не решить методами предыдущей главы, предназначенными для задач с начальными условиями. Далее для показа возможностей Mathcad будем использовать этот пример с R=1 и конкретным видом a(x)=const=1 и r(x)=const=0.1, описывающим случай изотропного (не зависящего от координаты х) рассеяния.  
 

ПРИМЕЧАНИЕ 1  

Модель, представленная на рис. 10.1, привела к краевой задаче для системы линейных ОДУ. Она имеет аналитическое решение в виде комбинации экспонент. Более сложные, нелинейные, задачи возможно решить только численно. Нетрудно сообразить, что модель станет нелинейной, если сделать коэффициенты ослабления и рассеяния зависящими от интенсивности излучения. Физически это будет соответствовать изменению оптических свойств среды под действием мощного излучения.  
 
 
 

ПРИМЕЧАНИЕ 2  

Модель встречных  световых пучков привела нас к  системе уравнений (10.1), в которые  входят производные только по одной  переменной х. Если бы мы стали рассматривать  более сложные эффекты рассеяния  в стороны (а не только вперед и  назад), то в уравнениях появились  бы частные производные по другим пространственным переменным у и z. В этом случае получилась бы краевая задача для уравнений в частных производных, решение которой во много раз сложнее ОДУ.

Для решения краевых  задач в Mathcad реализован наиболее популярный алгоритм, называемый методом стрельбы или пристрелки (shooting method). Он, по сути, сводит решение краевой задачи к решению серии задач Коши с различными начальными условиями. Рассмотрим сначала основной принцип алгоритма стрельбы на примере модели световых пучков (см. рис. 10.1), а встроенные функции, реализующие этот алгоритм, приведем в следующем разделе. 
 

Суть метода стрельбы заключается в пробном задании  недостающих граничных условий  на левой границе интервала и  решении затем полученной задачи Коши хорошо известными методами (см. главу 11). В нашем примере не хватает начального условия для y1(0), поэтому сначала зададим ему произвольное значение, например, у1(0)=10. Конечно, такой выбор не совсем случаен, поскольку из физических соображений ясно, что, во-первых, интенсивность излучения — величина заведомо положительная, и, во-вторых, отраженное излучение должно быть намного меньше падающего. Решение задачи Коши с помощью функции rkfixed приведено в листинге 10.1, причем в его последней строке выведены расчетные значения у0 и уг на правой крайней точке интервала х=1. Разумеется, они не совпадают, т. е. правое краевое условие не выполнено, из чего следует, что полученный результат не является решением поставленной краевой задачи.  

Листинг 10.1. Решение  пробной задачи Коши для модели (10.1)

   

График полученных решений показан на рис. 10.2 (слева). Из него также видно, что взятое наугад второе начальное условие не обеспечило выполнение граничного условия при  х=1. В целях лучшего выполнения этого граничного условия следует  взять большее значение y1(0), например, y1(0)=15, и вновь решить задачу Коши. Соответствующий результат показан  на том же рис. 10.2 (в центре). Граничное  условие выполняется с лучшей точностью, но опять-таки оказалось  недостаточным. Для еще одного значения y1(0)=20 получается решение, показанное на рис. 10.2 (справа). Из сравнения двух правых графиков легко заключить, что недостающее  начальное условие больше 15, но меньше 20. Продолжая подобным образом "пристрелку"  по недостающему начальному условию, возможно отыскать правильное решение краевой задачи.  

В этом и состоит  принцип алгоритма стрельбы. Выбирая  пробные начальные условия (проводя  пристрелку) и решая соответствующую  серию задач Коши, можно найти  то решение системы ОДУ, которое (с заданной точностью) удовлетворит граничному условию (или, в общем  случае, условиям) на другой границе  расчетного интервала.  

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

 

Рис. 10.2. Иллюстрация  метода стрельбы (продолжение листинга 10.1) 

Решение краевых  задач для систем ОДУ методом  стрельбы в Mathcad достигается применением двух встроенных функций. Первая предназначена для двухточечных задач с краевыми условиями, заданными на концах интервала:  

 sbval(z,x0,x1, D, load, score) — поиск вектора недостающих L начальных условий для двухточечной краевой задачи для системы N ОДУ:  
 

 z — вектор размера Lx1, присваивающий недостающим начальным условиям (на левой границе интервала) начальные значения;

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

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

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

 score (x1, у) — векторная функция размера Lx1, выражающая L правых граничных условий для векторной функции у в точке x1;

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

ВНИМАНИЕ!  

Решение краевых  задач в Mathcad, в отличие от большинства остальных операций, реализовано не совсем очевидным образом. В частности, помните, что число элементов векторов о и load равно количеству уравнений N, а векторов z, score и результата действия функции sbval— количеству правых граничных условий L. Соответственно, левых граничных условий в задаче должно быть (N-L) .  
 
 

Как видно, функция  sbval предназначена не для поиска собственно решения, т. е. неизвестных функций yi(x), а для определения недостающих начальных условий в первой точке интервала, т. е. yi(x0). Чтобы вычислить yi(х) на всем интервале, требуется дополнительно решить задачу Коши. Разберем особенности использования функции sbval на конкретном примере (листинг 10.2), описанном выше (см. разд. 10.1.1). Краевая задача состоит из системы двух уравнений (N=2), одного левого (L=1) и одного правого (N-L=2-1=1) граничного условия.  

Листинг 10.2. Решение  краевой задачи

 

Первые три строки листинга задают необходимые параметры  задачи и саму систему ОДУ. В четвертой  строке определяется вектор z. Поскольку правое граничное условие всего одно, то недостающее начальное условие тоже одно, соответственно, и вектор z имеет только один элемент z0. Ему необходимо присвоить начальное значение (мы приняли z0=10, как в листинге 10.1), чтобы запустить алгоритм стрельбы (см. разд. 10.1.2).  
 

ПРИМЕЧАНИЕ  

Начальное значение фактически является параметром численного метода и поэтому может решающим образом повлиять на решение краевой  задачи.  
 
 

В следующей строке листинга векторной функции load(x,z) присваиваются левые граничные условия. Эта функция аналогична векторной переменной, определяющей начальные условия для встроенных функций, решающих задачи Коши. Отличие заключается в записи недостающих условий. Вместо конкретных чисел на соответствующих местах пишутся имена искомых элементов вектора z. В нашем случае вместо второго начального условия стоит аргумент z0 функции load. Первый аргумент функции load — это точка, в которой ставится левое граничное условие. Ее конкретное значение определяется непосредственно в списке аргументов функции sbval. Следующая строка листинга определяет правое граничное условие, для введения которого используется функция score(х,у). Оно записывается точно так же, как система уравнений в функции о. Аргумент х функции score аналогичен функции load и нужен для тех случаев, когда граничное условие явно зависит от координаты х. Вектор score должен состоять из такого же числа элементов, что и вектор z.  

Реализованный в  функции sbval алгоритм стрельбы ищет недостающие начальные условия таким образом, чтобы решение полученной задачи Коши делало функцию score (х, у) как можно ближе к нулю. Как видно из листинга, результат применения sbval для интервала (0,1) присваивается векторной переменной и. Этот вектор похож на вектор z, только в нем содержатся искомые начальные условия вместо приближенных начальных значений, заданных в z. Вектор и содержит, как и z, всего один элемент и0. С его помощью можно определить решение краевой задачи у (х) (последняя строка листинга). Тем самым функция sbval сводит решение краевых задач к задачам Коши. График решения краевой задачи показан на рис. 10.3.  

На рис. 10.4 показано решение той же самой краевой  задачи, но с другим правым граничным  условием, соответствующим R=0, т. е. без  зеркала на правой границе. В этом случае слабый обратный пучок света  образуется исключительно за счет обратного  рассеяния излучения от лазера. Конечно, многие из читателей уже обратили внимание, что реальная физическая среда не может создавать такого большого рассеяния назад. Иными  словами, более реальны значения r(х)<<а (х). Однако когда коэффициенты в системе ОДУ при разных yi очень сильно (на порядки) различаются, система ОДУ становится жесткой, и функция sbval не может найти решения, выдавая вместо него сообщение об ошибке "Could not find a solution" (Невозможно найти решение).  
 

ВНИМАНИЕ!  

Метод стрельбы не годится  для решения жестких краевых  задач. Поэтому алгоритмы решения  жестких ОДУ в Mathcad приходится программировать самому (см. разд. 10.4). 

 

 

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

 bvaifit(z1,z2,x0,x1,xf,D,load1,load2,score) — поиск вектора недостающих граничных условий для краевой задачи с дополнительным условием в промежуточной точке для системы N ОДУ.  
 

z1 — вектор, присваивающий  недостающим начальным условиям  на левой границе интервала  начальные значения.

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