Кирьянов Д. - MathCad 11 (1077323), страница 51
Текст из файла (страница 51)
) - V.Yi= 001У^>2i+1L—»Д«^N-11XРис. 12.7. Сетка, покрывающая расчетный интервалПримечаниеСуществует множество способов аппроксимации дифференциальных уравнений разностными. От выбора конкретного варианта зависит не только простота,быстрота и удобство вычислений, но и сама возможность получения правильного ответа.Получилась система (по числу шагов) 2(ы-1) разностных линейных алгебраических уравнений с 2-N неизвестными Yi и yL. Для того чтобы она имелаединственное решение, надо дополнить число уравнений до 2-N. ЭТО МОЖНОсделать, записав в разностном виде оба граничных условия:Y0 = I 0 ,yN=R'YN.(2)Сформированная полная система алгебраических уравнений называетсяразностной схемой, аппроксимирующей исходную краевую задачу.
Обратите внимание, что правые части разностных уравнений системы (1) на каждом шаге записаны для левой границы шага. Такие разностные схемы называют явными, т. к. все значения Yi,! и y i + 1 находятся в левой частиуравнений. Полученную явную разностную схему легко записать в матричной формеA-z=B,(3)где z — неизвестный вектор, получающийся объединением векторов Y и у.Решив систему (3), мы получим решение краевой задачи.Глава 12. Краевые задачи313ПримечаниеНа самом деле, все несколько сложнее, поскольку, вообще говоря, необходимоеще доказать, что, во-первых, разностная схема действительно аппроксимирует дифференциальные уравнения и, во-вторых, при N ^ < » разностное решениедействительно сходится к дифференциальному.Процесс решения системы разностных уравнений называют также реализацией разностной схемы. Программа, которая решает рассматриваемую краевую задачу разностным методом, приведена в листинге 12.6.I Листинг 12.6.
Реализация явной разностной схемыа ( х ) : = 1г(х}:=0.01Ni := 1 ..i : = N + 1 .,AoANfi N-2-i := 1 ..1Ai , i + i := 1Ai , i + N := - Д - r ( i • Д - Д)2-N-li := - a ( i -Д - д ) - Д - 1=R := 1N-1A i , i := a ( i -Д - д ) - Д - 1AiЮ := 1 0 0IAi , i + i := 1Ai _ i_ N := Д • r ( i • Д - Д)1=1Ajg [j ;— —R2-N-lBi := 0Bo := 10A := s u b m a t r i x (A , 0 , 2 -N - 1 , 1 , 2 -N)U := l s o l v e (A, B)UM = 1 3 . 4 5 3UN_i = 3 1 .
7 8Дадим минимальные комментарии, надеясь, что заинтересовавшийся читатель с карандашом в руках разберется в порядке индексов и соответствииматричных элементов, а возможно, составит и более удачную программу.В первой строке листинга определяются функции и константы, входящие вмодель, во второй задается число точек сетки N=5 и ее равномерный шаг.Следующие две строки определяют матричные коэффициенты, аппроксимирующие уравнения для YL, а пятая и шестая — для у±.
Седьмая и восьмаястроки листинга задают, соответственно, левое и правое граничное условие,а строки с девятой по одиннадцатую — правые части системы (3). В следующей строке завершается построение матрицы А вырезанием из нее левого нулевого столбца. В предпоследней строке листинга применена встроенная функция lsolve для решения системы (3), а в последней выведенырассчитанные ею неизвестные граничные значения. Графики решения при-Часть HI. Численные методы314ведены на рис.
12.8, причем первые N элементов итогового вектора есть вычисленное излучение вперед, а последние N элементов — излучение назад.а.ч50К-1V1_**s»_—лл00<$•—0.5l1-ДРис. 12.8. Решение краевой задачи разностным методом(листинг 12.6)Как мы увидели, реализация в Mathcad разностных схем вполне возможна ине слишком трудоемка — предложенная программа состоит всего из двухдесятков математических выражений. Конечно, для их написания требуетсяи время, и часто кропотливые расчеты, но, собственно, в этом и состоитработа математика. Кстати говоря, при небольшом числе шагов, расчеты поразностным схемам не требуют существенного времени (программа, приведенная в листинге 12.6, работает быстрее, чем метод стрельбы, встроенныйв функцию sbvai). Существуют, кроме того, весьма очевидные для многихчитателей пути ускорения расчетов, связанные с применением более подходящих методов решения систем линейных уравнений с разреженной матрицей.12.3.2.
Жесткие краевые задачиОдин из случаев, когда применение разностных схем может быть очень полезным, связан с решением жестких краевых задач (подробнее о жесткихОДУ читайте в гл. 11). В частности, рассматриваемая задача о встречныхсветовых пучках становится жесткой при увеличении коэффициента ослабления а{х) в несколько десятков раз. Например, при попытке решить ее са(х) :=юо с помощью листинга 12.2 вместо ответа выдается сообщение обошибке "Can't converge to a solution.
Encountered too many integration steps"(He сходится к решению. Слишком много шагов интегрирования). Это инеудивительно, поскольку жесткие системы характерны тем, что требуютисключительно малого значения шага в стандартных алгоритмах.Для жестких задач неприменимы и явные разностные схемы, о которыхрассказывалось в предыдущем разделе (см. разд. 12.3.1). Результат расчетовГлава 12. Краевые задачи315по программе листинга 12.6, например с а(х) :=20 (рис. 12.9), дает характерную для неустойчивых разностных схем "разболтку" — колебания нарастающей амплитуды, не имеющие ничего общего с реальным решением.Рис. 12.9.
Неверное решение жесткой краевой задачипо неустойчивой явной разностной схемеВыходом из положения будет использование неявных разностных схем.Применительно к нашей задаче достаточно заменить правые части уравнений (I) значениями не на левой, а на правой границе каждого шага:ci+l> • Y i + 1K+ r< x i + l> ' У И(4)i + l) • Yi + i - r ( x i + l> ' Y i -Граничные условия, конечно, можно оставить в том же виде (2). Посколькумы имеем дело с линейными дифференциальными уравнениями, то и схему(4) легко будет записать в виде матричного равенства (3), перегруппировывая соответствующим образом выражение (4) и приводя подобные слагаемые. Разумеется, полученная матрица А будет иной, нежели матрица А ДЛЯявной схемы (1).
Поэтому и решение (реализация неявной схемы) можетотличаться от изображенного на рис. 12.9 результата расчетов по явной схеме. Программа, составленная для решения системы (4), приведена в листинге 12.7.Листинг 12.7. Реализация неявной разностной схемыдля жесткой краевой задачиа(х):=20N := 2 0i := 1 ..г{х)Д :=N - 1:-0.011N- 110 := 1 0 0R := 1•*. 1316AЧасть III. Численные методыi , i := 0 - 1i : = N + 1 ..Ai , i+i := a ( i -Д - д)-Д + 1Ai f i+H+l := -Д - r ( i -Д - д)2 -N-li , i := 0 - 1Ai (i+i := 1-а(л.-Д-д)-Дhi , i-N+i := ДAo , i : = 1AN, N•2i := 1 ..: =A^N, N: =~R2-N-lB i := ОB o := 1 0A := s u b m a t r i x ( A , 0 , 2 - N - 1 , 1 , 2 -N)U := l s o l v e ( A , B)UN=0.017U.-6N - 1 = 1 - 5 2 3 x 10He будем специально останавливаться на обсуждении листинга 12.7, поскольку он почти в точности повторяет предыдущий листинг.
Отличие заключается лишь в формировании матрицы А другим способом, согласнонеявной схеме. Решение, показанное на рис. 12.10, демонстрирует, что произошло "небольшое чудо": "разболтка" исчезла, а распределение интенсивностей стало физически предсказуемым. Обратите внимание, что (из-за взятого нами слишком большого коэффициента ослабления излучения) отраженный пучок света имеет очень маленькую интенсивность, и ее пришлосьпостроить на графике с увеличением в тысячу раз.i ;= 0 .N- 10 > + i -1000o--«>-oРис. 12.10. Решение краевой задачи разностным методом{листинг 12.7)ГЛАВА 13Дифференциальныеуравнения в частныхпроизводныхДифференциальные уравнения в частных производных представляют собой однуиз наиболее сложных и одновременно интересных задач вычислительной математики.
Эти уравнения характеризуются тем, что для их решения не существует единого универсального алгоритма и большинство задач требует своего собственного особого подхода. Уравнениями в частных производныхописывается множество разнообразных физических явлений, и с их помощью можно с успехом моделировать самые сложные явления и процессы(диффузия, гидродинамика, квантовая механика, экология и т. д.).Дифференциальные уравнения в частных производных требуют нахождения функции не одной, как для ОДУ, а нескольких переменных, напримерf(x,y) или f(x,t).
Постановка задач (см. разд. 13.1) включает в себя самоуравнение (или систему уравнений), содержащее производные неизвестнойфункции по различным переменным (частные производные), а также определенное количество краевых условий на границах расчетной области.Несмотря на то, что Mathcad обладает довольно ограниченными возможностями по отношению к уравнениям в частных производных, в нем имеетсянесколько встроенных функций, количество и возможности которых увеличились в новой 11-й версии (см.
разд. 13.3). Решать уравнения в частныхпроизводных можно и путем непосредственного программирования пользовательских алгоритмов (см. разд. 13.2). Автор совершенно сознательно сначала рассматривает численные методы для решения уравнений в частныхпроизводных, а уже затем описывает предназначенные для этого встроенные функции, чтобы читатель ясно осознавал, каким образом Mathcad 11производит расчеты. «Слепое» использование встроенных функций для решения уравнений в частных производных не всегда бывает успешным, и318Часть III.
Численные методыответственность за верный выбор их параметров часто ложится на пользователя, которому необходимо четко представлять основные принципы функционирования численных алгоритмов, примененных во встроенных функциях.Рассмотрение уравнений в частных производных построено на анализе различных вариаций уравнения диффузии тепла, которое допускает нагляднуюфизическую интерпретацию.13.1. Постановка задач13.1.1. Классификация уравненийв частных производныхПостановка задач для уравнений в частных производных включает определение самого уравнения (или системы нескольких уравнений), а также необходимого количества краевых условий (число и характер задания которыхопределяется спецификой уравнения).