Кирьянов Д. - MathCad 11 (1077323), страница 49
Текст из файла (страница 49)
гл. 11). В нашем примерене хватает начального условия для у х (0), поэтому сначала зададим емупроизвольное значение, например y ^ o j - i o . Конечно, такой выбор несовсем случаен, поскольку из физических соображений ясно, что, вопервых, интенсивность излучения — величина заведомо положительная, и,во-вторых, отраженное излучение должно быть намного меньше падающего.Решение задачи Коши с помощью функции rkfixed приведено в листинЛистинг 12.1.
Решение пробной задачи Коши для модели (12.1,1)г (х):= . 1а (х):= 1R :- 1- а (х}D ( х , у)-уо + г ( х ):=а(х)-ух - г ( х )уо\302Часть !П. Численные методыЮ := 100[1 0у :={ 10М := 10S •= r k f i x e d (у . 0 , 1 , М, D)График полученных решений показан на рис. 12.2 (слева). Из него видно,что взятое наугад второе начальное условие не обеспечило выполнение граничного условия при х=1. И понятно, что для лучшего выполнения этогограничного условия следует взять большее значение yi(0).
Возьмем, например, у!(0)=15, и вновь решим задачу Коши. Результат показан на том жерис. 12.2 (в центре). Граничное условие выполняется с лучшей точностью,но опять-таки оказалось недостаточным. Для еще одного значения у!(0)=20получается решение, показанное на рис. 12.2 (справа). Из сравнения двухправых графиков легко заключить, что недостающее начальное условиебольше 15, но меньше 20.
Продолжая подобным образом "пристрелку" понедостающему начальному условию, возможно отыскать правильное решение краевой задачи.В этом и состоит принцип алгоритма стрельбы. Выбирая пробные начальные условия (проводя пристрелку) и решая соответствующую серию задачКоши, можно найти то решение системы ОДУ, которое (с заданной точностью) удовлетворит фаничному условию (или, в общем случае, условиям)на другой фанице расчетного интервала.100(г)15.37100.510.51 010.5.
(о)Рис. 12.2. Иллюстрация метода стрельбы (листинг 12.1)Конечно, описанный алгоритм несложно запрограммировать самому, оформив его как решение системы заданных алгоритмически уравнений, выражающих граничные условия на второй границе, относительно неизвестныхпристрелочных начальных условий. Но делать этого нет необходимости, поскольку он оформлен в Mathcad в виде встроенных функций.Глава 12. Краевые задачи30312.1.3. Решение двухточечных краевых задачРешение краевых задач для систем ОДУ методом стрельбы в Mathcadдостигаетсяприменениемдвухвстроенныхфункций.Перваяпредназначена для двухточечных задач с краевыми условиями, заданнымина концах интервала.О sbvai(z,xO,xi,D, load, score) — поиск вектора недостающих L начальных условий для двухточечной краевой задачи для системы N ОДУ;•z — вектор размера LXI, присваивающий недостающим начальным условиям (на левой границе интервала) начальные значения;• хО — левая граница расчетного интервала;••xi — правая граница расчетного интервала;ioad(xO,z) — векторная функция размера NXI левых граничных условий, причем недостающие начальные условия именуются соответствующими компонентами векторного аргумента z;• score (xi, у) — векторная функция размера LXI, выражающая L правых граничных условий для векторной функции у в точке xi;• о{х,у) — векторная функция, описывающая систему ы ОДУ, размераNxi и двух аргументов — скалярного х и векторного у.
При этом у —это неизвестная векторная функция аргумента х того же размера NXI.Внимание!Решение краевых задач в Mathcad, в отличие от большинства остальных операций, реализовано не совсем очевидным образом. В частности, помните, чточисло элементов векторов D и l o a d равно количеству уравнений N, а векторовz, s c o r e и результата действия функции s b v a l — количеству правых граничных условий L.
Соответственно, левых граничных условий в задаче должно6blTb(N-L) .Как видно, функция sbval предназначена не для поиска собственно решения, т. е. неизвестных функций yi(x), а для определения недостающих начальных условий в первой точке интервала, т. е. yi(xO). Чтобы вычислитьyi(x) на всем интервале, требуется дополнительно решить задачу Коши.Разберем особенности использования функции sbval на конкретном примере (листинг 12.2), описанном выше (см. разд. 12.1.1). Краевая задача состоит из системы двух уравнений (N=2), ОДНОГО левого ( L = I ) И ОДНОГО правого ( N - L = 2 - I = I ) граничного условия.! Листинг 12.2.
Решение краевой задачиг ( х ) := . 1а ( х ) := 1R := 110 := 100''I304Часть III. Численные методыD ( х , у) :=( - а (х) ' у о + г (х){а (х) -у!-г (х)го := Юl o a d (хО , z) :=( 100zoI^)s c o r e ( x l , у) :- R 'Уо "~ У1I I := s b v a l ( z , 0 , 1 , D , l o a d , s c o r eI I = ( 18.555 >S :- r k f i x e d10, 1 , 10 , DПервые три строки листинга задают необходимые параметры задачи и самусистему ОДУ.
В четвертой строке определяется вектор z. Поскольку правоеграничное условие всего одно, то недостающее начальное условие тоже одно, соответственно, и вектор z имеет только один элемент z0. Ему необходимо присвоить начальное значение (мы приняли z o =io, как в листинге 12.1), чтобы запустить алгоритм стрельбы (см. разд 12.1.2).ПримечаниеНачальное значение фактически является параметром численного метода ипоэтому может сильно повлиять на решение краевой задачи.В следующей строке листинга векторной функции ioad{x,z) присваиваются левые граничные условия. Эта функция аналогична векторной переменной, определяющей начальные условия для встроенных функций, решающих задачи Коши.
Отличие заключается в записи недостающих условий.Вместо конкретных чисел на соответствующих местах пишутся имена искомых элементов вектора z. В нашем случае вместо второго начального условия стоит аргумент z0 функции load. Первый аргумент функции load — этоточка, в которой ставится левое граничное условие. Ее конкретное значениеопределяется непосредственно в списке аргументов функции sbval.Следующая строка листинга определяет правое граничное условие, для введения которого используется функция score{x,у). Оно записывается точнотак же, как система уравнений в функции D. Аргумент х функции scoreаналогичен функции load и нужен для тех случаев, когда граничное условиеявно зависит от координаты х. Вектор score должен состоять из такого жечисла элементов, что и вектор z.Реализованный в функции sbval алгоритм стрельбы ищет недостающие начальные условия таким образом, чтобы решение полученной задачи Кошиделало функцию score (х,у) как можно ближе к нулю.
Как видно из листинга, результат применения sbval для интервала ( о , и присваивается век-Глава 12. Краевые задачи305торной переменной и. Этот вектор похож на вектор z, только в нем содержатся искомые начальные условия вместо приближенных начальных значений, заданных в z. Вектор и содержит, как и z, всего один элемент и 0 .С его помощью можно определить решение краевой задачи у(х) (последняястрока листинга). Тем самым, функция sbval сводит решение краевых задачк задачам Коши. График решения краевой задачи показан на рис.
12.3,Рис. 12.3. Решение краевой задачидля R=l (листинг 12.2)Рис. 12.4. Решение краевой задачидля R=0На рис. 12.4 показано решение той же самой краевой задачи, но с другимправым граничным условием, соответствующим R=O, Т. е. без зеркала направой границе. В этом случае слабый обратный пучок света образуется исключительно за счет обратного рассеяния излучения от лазера. Конечно,многие из читателей уже обратили внимание, что реальная физическая среда не может создавать такого большого рассеяния назад.
Иными словами,более реальны значения г ( х ) « а ( х ) . Однако, когда коэффициенты в системе ОДУ при разных у± очень сильно (на порядки) различаются, системаОДУ становится жесткой, и функция sbval не может найти решения, выдавая вместо него сообщение об ошибке("СоиМ not find a solution").Внимание!Метод стрельбы не годится для решения жестких краевых задач. Поэтому алгоритмы решения жестких ОДУ в Mathcad приходится программировать самому(см. разд. 12.3).12.1.4. Решение краевых задач с дополнительнымусловием в промежуточной точкеИногда дифференциальные уравнения определяются с граничными условиями не только на концах интервала, но и с дополнительным условием в306Часть III.
Численные методынекоторой промежуточной точке расчетного интервала. Чаще всего такиезадачи содержат данные о негладких в некоторой внутренней точке интервала решениях. Для них имеется встроенная функция bvaifit, также реализующая алгоритм стрельбы.Оbvaifit(zl,z2,xO,xl,xf,D,loadl,Ioad2,score) —ПОИСКвекторане-достающих граничных условий для краевой задачи с дополнительнымусловием в промежуточной точке для системы N ОДУ;• zi — вектор, присваивающий недостающим начальным условиям налевой границе интервала начальные значения;• z2 — вектор того же размера, присваивающий недостающим начальным условиям на правой границе интервала начальные значения;• хо — левая граница расчетного интервала;• xi — правая граница расчетного интервала;•xf — точка внутри интервала;• о(х,у) — векторная функция, описывающая систему N ОДУ, размераNxi и двух аргументов — скалярного х и векторного у.
При этом у —это неизвестная векторная функция аргумента х того же размера NXI;•loadi(xO,z) — векторная функция размера NXI левых граничных условий, причем недостающие начальные условия поименовываютсясоответствующими компонентами векторного аргумента z;• ioad2(xi,z) — векторная функция размера NXI правых граничных условий, причем недостающие начальные условия поименовываютсясоответствующими компонентами векторного аргумента г;•score(xf,y) — векторная функция размера NXI, выражающая внутреннее условие для векторной функции у в точке xf.Обычно функция b v a i f i t применяется для задач, в которых производнаярешения имеет разрыв во внутренней точке xf. Некоторые из таких задачне удается решить обычным методом пристрелки, поэтому приходится вестипристрелку сразу из обеих граничных точек.
В этом случае дополнительноевнутреннее условие в точке xf является просто условием сшивки в ней левого и правого решений. Поэтому для данной постановки задачи оно записывается в форме scorefxf, у) :=у.Рассмотрим действие функции bvaifit на знакомом примере модели взаимодействия пучков света (см. рис. 12.1), предположив, что в промежуткемежду xf=o.5 и xi=i.O находится другая, оптически более плотная среда сдругим коэффициентом ослабления излучения а (х) =з.
Соответствующаякраевая задача решена в листинге 12.3, причем разрывный показатель ослабления определяется в его второй строке.Глава 12. Краевые задачи307Листинг 12.3. Краевая задача страничным условиемв промежуточной точке10 := 100ах:=г (х) := .11R := 1if х < 0.53 otherwiseD(x, у) :=-а (х) -уо + г (х) -У1а (х)- г (х) -уоzl 0 := 20100loadl (xO , zl) :=z2 0 := 30Ioad2 (xl , z2) :=(z2 0R-z2,score ( xf , у) := уII :=bvalfit ( zl , z2 , 0 , 1, 0. 5 , D, loadl , Ioad2 , score)5.618II =13.801Система уравнений и левое краевое условие вводится так же, как и в предыдущем листинге для функции sbvai. Обратите внимание, что таким жеобразом записано и правое краевое условие.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
