Кирьянов Д. - MathCad 11 (1077323), страница 44
Текст из файла (страница 44)
п.Поскольку решение вторым способом одного ОДУ мало чем отличается отрешения систем ОДУ (см. разд. 11.3), приведем пример его использованияв задаче из листинга 11.1 практически без комментариев (см. листинг 11.2)и с помощью одной из трех существующих для этих целей встроенныхфункций rkfixed.
Обратите внимание только на необходимость явного задания количества точек интегрирования ОДУ м=юо в третьей строке листинга, а также на получение результата, в отличие от вычислительного блока, не в виде функции, а в виде матрицы размерности мх2. Она состоит издвух столбцов: в одном находятся значения аргумента t (от to до t i включительно), а в другом соответствующие значения искомой функции у ( t ) .Листинг 11.2. Решение задачи Коши для ОДУ первого порядкавторым способому := 0 . 1D ( t , у) : = у - у ^м := 100у := r k f i x e d ( y , 0 , 1 0 , M , D )СоветВ листинге 11.2 приведен пример не лучшего стиля Mathcad-программирования. Сначала переменной у присвоено значение скаляра у=0.1, а затем этойже переменной присвоено матричное значение (результат решения ОДУ).
Старайтесь избегать такого стиля, который ухудшает читаемость программы и может приводить, в более сложных случаях, к трудно опознаваемым ошибкам.Неплохим решением было бы назвать результат по-другому, например и.Глава 11, Обыкновенные дифференциальные уравнения271График решения рассматриваемого уравнения показан на рис. 11.1. Обратите внимание, что он соответствует получению решения в матричном виде(листинг 11.2), поэтому по осям отложены соответствующие столбцы, выделенные из матрицы у оператором <>.ПримечаниеПример, решенный в листингах 11.1—11.2, взят из области математическойэкологии и описывает динамику популяций с внутривидовой конкуренцией.Сначала происходит рост численности популяции, близкий к экспоненциальному, а затем выход на стационарное состояние.Рис.
11.1. Решение уравнения у ' =у-у 2 (листинг 11.2)11.2. ОДУ высшего порядкаОбыкновенное дифференциальное уравнение с неизвестной функцией y ( t ) ,(К)в которое входят производные этой функции вплоть до у (t), называетсяОДУ N-ГО порядка. Если имеется такое уравнение, то для корректной постановки задачи Коши требуется задать ы начальных условий на саму функцию y(t) и ее производные от первого до (и-и-го порядка включительно.В Mathcad 11 можно решать ОДУ высших порядков как с помощью вычислительного блока Given/odesolve, так и путем сведения их к системамуравнений первого порядка.Внутри вычислительного блока:П ОДУ должно быть линейно относительно старшей производной, т.
е.фактически должно быть поставлено в стандартной форме;•начальные условия должны иметь форму y(t)=b или y w (t)=b, а не более сложную (как, например, встречающаяся в некоторых математических приложениях форма у (t) +у' (t) =ь).Часть III. Численные методы272В остальном, решение ОДУ высших порядков ничем не отличается от решения уравнений первого порядка (см. разд. 11.1), что иллюстрируется листингом 11.3. Как Вы помните, допустимо написание производной как в виде знака дифференциала (так в листинге 11.3 введено само уравнение), таки с помощью штриха (так введено начальное условие для первой производной).
Не забывайте пользоваться булевыми операторами при вводе уравнений и начальных условий. Полученное решение y(t) показано на рис. 11.2.Листинг 11.3. Решение задачи Коши для ОДУ второго порядкаGivenЛy(t)dtу ( 0 ) = 0.1- — y(t)a tу ' (0) = 0у := Odesolve ( t , 50Рис. 1 1 .
2 . Решение уравнения осциллятора (листинг 11.3)ПримечаниеВ листинге 11.3 решено уравнение затухающего гармонического осциллятора,которое описывает, например, колебания маятника. Для модели маятника y ( t )описывает изменения угла его отклонения от вертикали, y ' ( t ) — угловую скорость маятника, y " ( t ) —ускорение, а начальные условия, соответственно, начальное отклонение маятника у (0) =0.1 и начальную скорость у' (0) = 0.Второй способ решения ОДУ высшего порядка связан со сведением его кэквивалентной системе ОДУ первого порядка.
Покажем на том же примереиз листинга 11.3, как это делается. Действительно, если формально обозна-Глава 11. Обыкновенные дифференциальные уравнения273чить y o ( t ) = y ( t ) , a y i ( t ) s y ' ( t ) = y 0 ' ( t ) , т о исходное у р а в н е н и е з а п и ш е т с ячерез ф у н к ц и и у 0 ( t ) и у г ( t ) в виде системы двух ОДУ:|Уо"=У1.( 1 )[ y i ' + O . l • Yi + 1 • У 0 = 0.Именно эта система решается в качестве примера в разд. 11.3. Таким образом, любое ОДУ N-ГО порядка, линейное относительно высшей производной, можно свести к эквивалентной системе N дифференциальных уравнений.11.3. Системы ОДУ первого порядкаMathcad требует, чтобы система дифференциальных уравнений была представлена в стандартной форме:Уо1( t ) = fo( yyN-i(t) , t ) ,( t ) , y i ( t ) , . .
. ,y( t ) ,Yi(t) ,H-i(t) ,t) ,0Yi' ( t ) = f i ( yoУ»-!1.(О{ t ) = fM1( y0( t }, y i ( t },yH-i(t) , t )Задание системы (1) эквивалентно следующему векторному представлению:У <t)=F<Y(t> , t ) ,(2)где Y и У — соответствующие неизвестные векторные функции переменнойt размера NXI, а р — векторная функция того же размера и {N+1} количества переменных (N компонент вектора и, возможно, t).
Именно векторноепредставление (2) используется для ввода системы ОДУ в среде Mathcad.Для того чтобы определить задачу Коши для системы ОДУ, следует определить еще N начальных условий, задающих значение каждой из функцийYi(tO) в начальной точке интегрирования системы to. В векторной формеони могут быть записаны в видеY(tO)=B,(3)где в — вектор начальных условий размера NXI, составленный изПримечаниеКак Вы заметили, задача сформулирована для систем ОДУ первого порядка.Однако если в систему входят и уравнения высших порядков, то ее можно свести к системе большего числа уравнений первого порядка, подобно тому как этобыло сделано на примере уравнения осциллятора (см.
разд. 11.2).Обратите внимание на необходимость векторной записи как самого уравнения, так и начального условия. В случае одного ОДУ соответствующие векторы имеют только один элемент, а в случае системы N>I уравнений — N.274Часть III. Численные методы11.3.1. Встроенные функциидля решения систем ОДУВ Mathcad 11 имеются три встроенные функции, которые позволяют решатьпоставленную в форме (2—3) задачу Коши различными численными методами.Пrkf ixedtyO, t o , ti,M,D) — метод Рунге-Кутты с фиксированным шагом;ПRkadapt (yo, t o , ti,M,D) — метод Рунге-Кутты с переменным шагом;•B u i s t o e r ( у - о , t o , t i , M , D ) — метод Булирша-Штера;•уо — вектор начальных значений в точке to размера NXI;•to — начальная точка расчета;•t i — конечная точка расчета;•м — число шагов, на которых численный метод находит решение;•D — векторная функция размера NXI двух аргументов — скалярного tи векторного у.
При этом у — искомая векторная функция аргументаt того же размера NXI.Внимание!Соблюдайте регистр первой буквы рассматриваемых функций, поскольку этовлияет на выбор алгоритма счета, в отличие от многих других встроенныхфункций Mathcad, например Find=find (см. разд. 11.3.2).Каждая из приведенных функций выдает решение в виде матрицы размера(M+i)x(N+i).
В ее левом столбце находятся значения аргумента t, делящиеинтервал на равномерные шаги, а в остальных N столбцах — значения искомых функций y o ( t ) ,yl{t), . . . , y N - i ( t ) , рассчитанные для этих значенийаргумента. Поскольку всего точек (помимо начальной) м, то строк в матрице решения будет всего м+1.В подавляющем большинстве случаев достаточно использовать первуюфункцию rkfixed, как это показано и листинге 11.4 на примере решениясистемы ОДУ модели осциллятора с затуханием (см.
разд. 11.2).Листинг 11.4. Решение системы двух ОДУD ( t , у ) :=-УО - 0 . 1уО :=М := 1 0 0и : - r k f i x e d ( у О , О , 5 0 , М , D)Глава 11. Обыкновенные дифференциальные уравнения275Самая важная — это первая строка листинга, в которой, собственно, определяется система ОДУ. Сравните рассматриваемую систему (разд. 11.2,1),записанную в стандартной форме, с формальной ее записью в Mathcad, чтобы не делать впоследствии ошибок. Во-первых, функция D, входящая вчисло параметров встроенных функций для решения ОДУ, должна бытьфункцией обязательно двух аргументов. Во-вторых, второй ее аргумент должен быть вектором того же размера, что и сама функция в. В-третьих, точнотакой же размер должен быть и у вектора начальных значений уо (онопределен во второй строке листинга).Не забывайте, что векторную функциюD(t,y) следует определять через компоненты вектора у с помощью кнопкинижнего индекса (Subscript) с наборнойпанели Calculator (Калькулятор) или нажатием клавиши <[>.
В третьей строкелистинга определено число шагов, накоторых рассчитывается решение, а егопоследняя строка присваивает матричнойпеременной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (о, 50).:•0- I"'•\.1 •••"••г- А0 01 G05 ООВВ -0 0472! 0 056 -0 08а 15 ООП -0 093•4 2 -0 033 -0 0825' 25 -0 06G -0 053% 3 -G0B4 -0 013" -I 7 35 -0 0В 0 029 ,:8•1-0 057 0 062.4' •15 -0 021 0 078W5 001В 0 07511 55 0 051 0 0541'/60 021га 65 ЮМ0 071^015и7 0 056-0 04S35 75 0D2B•0 064•у•Как выглядит все решение, показано нарис.
11.3. Размер полученной матрицыбудет равен ( M + D X ( N + D , Т. е. ioix3.Рис. 1 1 . 3 . Матрица решенийПросмотреть все компоненты матрицы и,системы уравнений (листинг 11.4)которые не помещаются на экране, можно с помощью вертикальной полосыпрокрутки. Как нетрудно сообразить, на этом рисунке отмечено выделениемрасчетное значение первого искомого вектора у0 на 12-м шаге и 12(1 =о.О7.Это соответствует, с математической точки зрения, найденному значениюУо (б .0) =0 .07. Для вывода элементов решения в последней точке интервалаиспользуйте выражения типа иМ|1=7 .523хю~3.Внимание!Обратите внимание на некоторое разночтение в обозначении индексов вектораначальных условий и матрицы решения.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
