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

Рис. 10.8. Решение краевой  задачи разностным методом (продолжение  листинга 10.6)

 
 

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

Один из случаев, когда применение разностных схем может  быть очень полезным, связан с решением жестких краевых задач (подробнее  о жестких ОДУ читайте в  главе 11). В частности, рассматриваемая задача о встречных световых пучках становится жесткой при увеличении коэффициента ослабления а(х) в несколько десятков раз. Например, при попытке решить ее с а(х) :=100 с помощью листинга 10.2 вместо ответа выдается сообщение об ошибке "Can't converge to a solution. Encountered too many integration steps" (He сходится к решению. Слишком много шагов интегрирования). Это и неудивительно, поскольку жесткие системы характерны тем, что требуют исключительно малого значения шага в стандартных алгоритмах.  

Для жестких задач  неприменимы и явные разностные схемы, о которых рассказывалось в предыдущем разделе (см. разд. 10.3.1). Результат расчетов по программе листинга 10.6, например с а(х) :=20 (рис. 10.9), дает характерную для неустойчивых разностных схем "разболтку" — колебания нарастающей амплитуды, не имеющие ничего общего с реальным решением.

 
 

Рис. 10.9. Неверное решение  жесткой краевой задачи по неустойчивой явной разностной схеме  
 

Выходом из положения  будет использование неявных  разностных схем. Применительно к  нашей задаче достаточно заменить правые части уравнений (10.1) значениями не на левой, а на правой границе каждого  шага:

   

Граничные условия, конечно, можно оставить в том  же виде (10.5). Поскольку мы имеем дело с линейными дифференциальными  уравнениями, то и схему (10.7) легко  будет записать в виде матричного равенства (10.6), перегруппировывая соответствующим  образом выражение (10.7) и приводя  подобные слагаемые. Разумеется, полученная матрица А будет иной, нежели матрица А для явной схемы (10.4). Поэтому и решение (реализация неявной схемы) может отличаться от изображенного на рис. 10.9 результата расчетов по явной схеме. Программа, составленная для решения системы (10.7), приведена в листинге 10.7.  

Листинг 10.7. Реализация неявной разностной схемы для  жесткой краевой задачи

   

Не будем специально останавливаться на обсуждении листинга 10.7, поскольку он почти в точности повторяет предыдущий листинг. Отличие  заключается лишь в формировании матрицы А другим способом, согласно неявной схеме. Решение, показанное на рис. 10.10, демонстрирует, что произошло "небольшое чудо": "разболтка" исчезла, а распределение интенсивностей стало физически предсказуемым. Обратите внимание, что (из-за взятого нами слишком большого коэффициента ослабления излучения) отраженный пучок света имеет очень маленькую интенсивность, и ее пришлось построить на графике с увеличением в тысячу раз.  
 

 

Рис. 10.10. Решение  краевой задачи по неявной разностной схеме (продолжение листинга 10.7) 

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

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

a(x,Y,y)=a0(x)-ε(x)(Y+y), (10.8)  

где а0(х) — постоянная (не зависящая от интенсивности света) составляющая коэффициента поглощения, а е(х) — малый коэффициент пропорциональности. В этом случае модель (10.1) перепишется в виде:

   

Таким образом, физическая подоплека квадратичной нелинейности задачи (10.9) связана с зависимостью коэффициента поглощения среды от Y+y. Иными словами, чем мощнее свет в некоторой точке среды, тем меньше локальный коэффициент поглощения среды в этой точке.  

Механизм этого  явления будет хорошо понятен, если предположить, что средой является воздух, насыщенный парами жидкости (попросту говоря, туман). Чем сильнее в некоторой  точке туман, тем больше локальное  поглощение. Если свет, проходящий через  туман, имеет малую интенсивность, то он не оказывает на среду никакого (или почти никакого) воздействия. Однако если пучок света будет  очень мощным, за счет его поглощения пары воды будут разогреваться и, следовательно, испаряться. По мере испарения  капель воды туман станет рассеиваться, и среда просветлеет (т. е. ее коэффициент  поглощения уменьшится). Таким образом, чем мощнее свет, тем активнее пойдет процесс испарения, тем существеннее уменьшится коэффициент поглощения. Именно такой эффект, совершенно понятный с физической точки зрения, и описывает  формула (10.8) и, соответственно, модель (10.9).  

Подчеркнем еще  раз и физический смысл коэффициентов  а0(х) и ε(х). Первый описывает постоянный (не зависящий от интенсивности света) коэффициент поглощения среды, а второй — его зависимость от эффекта разогрева среды светом. Чем больше коэффициент а0(х), тем более жесткой является краевая задача; а чем больше ε(х), тем сильнее ее нелинейность. В предельном случае ε(х)->0 задача (10.9) переходит в задачу (10.1), становясь линейной. 

Встроенная функция  sbval, реализующая в Mathcad алгоритм стрельбы (см. разд. 10.2), может справляться и с нелинейными задачами. Приведем пример решения краевой задачи (10.9) (листинг 10.8 и рис. 10.11) с теми же граничными условиями, что были поставлены для модели (10.1). Интерес представляет решение у, описывающее рост интенсивности отраженного пучка по мере его распространения справа-налево (обратите внимание, что на рис. 10.11 оно отложено с десятикратным увеличением).

Листинг 10.8. Решение  нелинейной краевой задачи методом  стрельбы

 

Следует помнить, что  чем сильнее нелинейность и чем  жестче краевая задача, тем более  сильные требования предъявляются  к начальному значению алгоритма, ввод которого осуществляется посредством  функции load. Попробуйте повторить расчеты листинга 10.8 с иным сочетанием параметров а0(х) и ε(х), и вы увидите, что с немного более жесткими задачами алгоритм стрельбы уже не справляется.  
 

ПРИМЕЧАНИЕ  

На компакт-диске, прилагаемом к книге, вы найдете  еще более интересное решение  для другой задачи, с нелинейностью  типа (y+Y)2.  

 
 
 

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

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

Подробная разработка алгоритмов и соответствующих программ Mathcad выходит далеко за пределы данной книги, поэтому мы конспективно представим один из приемов решения нелинейных краевых задач, сводящийся к их линеаризации. В общих чертах подход заключается в следующем. Предположим, что некоторое приближение (обозначим его Y(х) и у(х)) к решению нелинейной задачи (10.9) нам известно, и можно считать, что Y->Y+Z и y->y+z, где z и z — близкие к нулю функции х. Тогда, пользуясь их малостью (по сравнению с Y0 и у0), можно разложить нелинейные слагаемые в уравнениях (10.9) в ряд Тейлора по z и z. Получим:  

Y'+Z' = -aY + ry + ε(Y2 + Yy) - aZ + rz + 2εY(Z + z);     (10.10) 

y'+z' = ay - rY - ε(y2 + Yy) + az - rz - 2εy(Z + z) .  

Теперь, поскольку Y(x) и у(х) являются приближением к решению исходной задачи, то можно считать, что они (приблизительно) удовлетворяют и уравнению, и граничным условиям. Тогда, вычитая (10.9) из (10.10), получим краевую задачу для z (х) и z (х):  

Z' = -aZ + rz + 2εY(Z + z);  

z' = az - rZ - 2εy(Z + z) ; (10.11)  

Z(0) = z(0) = Z(l) = z(l) = 0.  

Это и есть та самая  новая задача для "маленьких" функций  z (х) и z (х), которую надо решить, и которая, благодаря малости z и z, является линейной. Вся беда в том, что мы не знаем "больших" функций Y(x) и у(х), а они, увы, входят в задачу (10.11). Тем не менее рецепт получения этих функций напрашивается сам собой: если нелинейность исходных ОДУ не слишком сильная, то в качестве Y и у можно взять решение линейной краевой задачи, т. е. задачи (10.9) с ε(х)=0 (см. разд. 10.4.1).  

Сказанное иллюстрируют листинги 10.9 и 10.10, первый из которых  решает линейную краевую задачу (являясь, фактически, небольшой модификацией листинга 10.7 из разд. 10.4.2), а второй решает линеаризованную задачу (10.11), учитывая результат листинга 10.9. В листинге 10.9 матрица о и вектор правых частей в являются разностной аппроксимацией ОДУ (ее первые N строк аппроксимируют первое уравнение, а оставшиеся N строк — второе). Такой же смысл и точно такую же структуру имеют матрица с и вектор правых частей с для второй (линеаризованной) задачи (10.11). Для решения сеточных уравнений Dу=B и Cz=G используется (конечно, весьма неэкономично) встроенная функция isolve, реализующая алгоритм Гаусса.  

Важно привлечь внимание читателя к последним строкам  листинга 10.9. В них осуществляется интерполяция полученного решения  системы сеточных уравнений для  того, чтобы в нелинейной задаче (в листинге 10.10) можно было использовать непрерывные "большие" функции  из линейной задачи. В последней  строке листинга 10.10 осуществляется сложение "больших" и "маленьких" функций (результатов листингов 10.9 и 10.10) для  получения полного решения нелинейной задачи (10.9), которое изображено на рис. 10.12.  

Не будем давать дополнительных комментариев к Mathcad-программам, надеясь, что читатель, заинтересовавшийся нелинейным примером со световыми пучками  и эффектом разогрева светом среды, сам разберется в листингах, тем  более что техника разработки разностных схем была нами детально разобрана  раньше (см. разд. 10.4).  

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

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

 

 

Последний важный момент, который следует обозначить, связан с решением задач, обладающих значительной нелинейностью. Решение, приведенное  в листингах 10.9, 10.10 и на рис. 10.12, согласно самой постановке, должно не сильно отличаться от решения линейной краевой  задачи, поскольку функции z(x) и z(x) малы по сравнению с Y(x) и у(х). Если же нелинейность сильная, то решение может сильно отличаться от Y и у, и линеаризация (10.11) будет просто неправильной. В этом случае следует слегка усложнить алгоритм решения нелинейной краевой задачи.  

 

Рис. 10.12. Решение  нелинейной краевой задачи разностным методом (продолжение листингов 10.9 и 10.10)  
 

Обозначим полученное в результате решение, как и в  листинге 10.10, вектором J(ε), подчеркивая тем самым его зависимость от параметра нелинейности. Очевидно, что J(0) есть решение линейной задачи. Для того чтобы решить задачу с сильной нелинейностью, т. е. довольно большим ε=ε1, можно организовать продолжение по ε как по параметру. Иными словами, используя в качестве начального приближения J (0), можно решить задачу для другого, малого ε=Дε, получив J(Aε), затем, взяв это J(ε) в качестве приближенного решения, получить J(2Aε) и т. д. малыми шагами добраться до желаемого ε1.  
 

ПРИМЕЧАНИЕ  

Упрощенную реализацию этого алгоритма вы найдете на компакт-диске, прилагаемом к книге. Она связана с выводом во внешний  файл данных результата задачи из листинга 10.10 и считыванием из него же этих данных в качестве входной информации для следующей итерации. В качестве нулевой итерации используется решение  линейной задачи, выводимое предварительно в файл из усовершенствованного листинга 10.9.