Говорухин В., Цибулин Б. Компьютер в математическом исследовании (1185927), страница 49
Текст из файла (страница 49)
В настоящее время развита соответствующая теория и известны многочисленные реализации методов, см. ссылки в (71]. Разберем непосредственный вывод исходя из общей схемы метода и ограничимся получением формулы для двухстадийного метода. Отметим, что в программе для краткости отсутствует какая-либо оптимизация при получении и последующем решении системы уравнений для искомых коэффициентов. Будем приводить команды Мар!е в математической и обычной нотациях.
Напомним, что в Мар!е для преобразования уже набранных строк ввода к нужной нотации достаточно нажать (отжать) значок Х на панели Соптехт Ваг. В математической нотации оператор суммирования записывается в виде суммы, а индексы обретают естественный вид.
260 Глава 1О. Примеры решения задач В начале программы полезно поставить команду гезтагс, чтобы очистить все назначения. Определим переменную п для указания числа этапов (стадий): > гехгагнл» 2 Подготовим формулу для вычисления решения в точке ь-Ь: > Арр:ю у + ( Х 'Ы5' () Арр:=у+ Ь, Ь, +Ь,Ь, Для сравнения покажем, как данная команда выглядит в обычной нотации: > Арр."у>ЗОМ("Ь(1)*И(1]".3-1..П); Вспомогательные функции для вычисления стадий поместим в специальную переменную ии (последовательность выражений): ю — ! > хх:= Ь Г Л у е ~~~ 'а х' у = [вел(л — т,т = О ..л — 1)] 1 ЬЬ:=Рл=ЬГ(у+а,,я,),Ь, =ЬГ(у) При п-2 число неизвестных коэффициентов равно трем; > чаг ч> зсд(Ь„У = 1 .. л ), зсд(зсд(а, „! = 1,.У вЂ” ! )У = 1 ., л ) з чаг; Ь, Ь, а Запишем исходное дифференциальное уравнение: > Ие чч — у(х) = Г(у(х)) д дх Сформируем последовательность выражений для производных высоких порядков функции у(х), полученных в предположении, что у(х) является решением данного дифференциального уравнения.
Для этого в цикле дифференцируем уравнение и используем соотношения, накапливаемые в переменной вез, для преобразования производных в правой части: > дез:-Ое: гог 2 !о п оо Ое: -Ш Гц)пз(ое), х)-з воз (дев, О1 ( (( гпз Ые) . х) ): оез:-оез,се: об: Вывод сформированной последовательности производных оез опустим. Например, последний элемент переменной оез дает выражение третьей производной функции у(х) через производные правой части уравнения ое: > вез[-П: дз дх' д „„,, ",„,х,),„х),, у(„)),у(„» Для точного решения найдем разложение в ряд по Ь: > Ехасг вч сопчса(яснев(у(х+ Ь ), Ь.
л + 2), д((Г) Вывод формуя явного метода Рунге-Купы 261 Е л.ау(,)+~,(„)~И+ Г,(,»И2, ~,(„)'Из+0(И ) Используя подготовленную переменную Оез, произведем замену производных: > ЕО;= виЬя(у(х) =у, тиьв((Мех), Ехасг)) ЕО:= у+((У)И+ Р(Г)(у) Г(у) И + -(Р )(()(У) Ку) +-Р(Г)(у) Г(у) )И +0(И ) 2 (6 б Теперь разложим в ряд по Ь приближенное решение в точке Г+Ь: > АО;= яепее(яиьа(Ыс,Арр), И, а+2) АО х У+ ( Ь, ((У) + Ь, Г(у) ) И+ Ь, Р(гКУ) ат, ((У) И'+- Ь, (Р )(Г)(у) а,, ((У)' Ит + 0(И4) Вычислим невязку и отбросим малые величины порядка и+2, превратив выражение в полипом: > иг;= ЯнпР(КУ(соптеп(АΠ— еО, ро(унон)) В дальнейшем удобно рассматривать производные разных порядков как некоторые скалярные величины.
Для этого подготовим ряд вспомогательных замен: > аиЬ.-яо)((099))( 0 (у)-Г(2).)-з..п); хаЬ:= ((у) = Р, Р(г')(у) = Р, (Р )(г)(у) = Р Запишем последовательность введенных вспомогательных величин: > со Е: "йдр( гп5. (аиЬ) ) П ' со Р:а Ра Рг Е Используя эти обозначения, перепишем выражение невязки, сгруппировав коэф- фициенты при степенях Ь и вспомогательных величинах Г,, Е,, Гк > и(0:= сонес((яиья(хиь, (Р ).
(И, со р), йхгг(Ьи(еО ) ч(О:= -""" 1"" ""'( "'" 1"""'"-' "'"" Из полученного выражения для невязки видно, что никаким подбором коэффициентов не удастся обнулить второе слагаемое. Таким образом, для двухстадийного метода максимум возможного — это.метод второго порядка точности, который получится, если найдутся значения коэффициентов, при которых пропадут члены первого и второго порядка по Ь.
Соберем коэффициенты при различных степенях шага Ь и учтем, что при одной степени Ь могут быть различные комбинации производных (у нас это теперь выражения, составленные из элементов со Р). Выделение уравнений помогает осуществить следующий фрагмент программы, где используется команда суре для определения типа выражения и превращения суммы в последовательность выражений: 262 Глава 10. Примеры решения задач > ео: Нцма Гог З Со и оо р:-сое(((рвго,ь"З); !( Суре(р,".") Снеп тог в гп р Оо ее:-ее.ор(1,в): ое: е15е ее:-ее.ор(1.р); (1; оо: В результате уравнения будут собраны как элементы переменной ец.
Выведем число уравнений и их вид: > поря([еа)); еа 2 1 Ь вЂ” 1+Ь,— — +Ь а Остается решить полученную систему уравнений относительно неизвестных ко- эффициентов: > ю:а во!че( ( еа ), ( >аг ) ) 1 1 ю;=(а =- —,Ь =Ь,Ь 1 — Ь ) г.! 2Ь' г г'! г Подставим найденные величины в невязку: > в!пгр!!Гу(воЬв(еа, г(г)) 1 Ь'((У) (-3 (() )(Г)(У) ((У) + 4 (() )(Г)(У) ((У) Ьг+ 4 1)(Г)(У)' Ьг) 24 Ь Итак, невязка имеет третий порядок по Ь.
Результирующие формулы получим подстановкой. Вспомогательные функции после подстановки коэффициентов примут вид (пустые скобки позволяют превратить список в последовательность выраженийу: БпЬв(еа, (ае]) 1! 1 2 г У ' ! г Таким образом, для вычисления значения функции у(С) за шаг имеем однопара- метрическое семейство формул: » впЬв(ю,Арр) У+(-Ь,+ !) Я, +Ь,а, Итак, получено семейство формул Рунге-Кутты второго порядка с одним свободным параметром. Всяи рассмотреть формулы Рунге-Кутты с четырьмя стадиями (и 4) „,то для метода четвертого порядка точности (выкладки опускаемя воаучается систевеа де)йи Вывод формул явного негода Рун?е-Купа! '2бмв уравнений с десятью неизвестными параметрами.
Известно множество вариантов данного метода, в частности классический метод Рунге-Кутты и так называемый метод «три восьмых» 171]. Анализируя полученную систему уравнений, можно указать еще одно решение следующего вида: 1 3 3 1 1 -1 а =-,Ь "--,Ь =-,Ь =-,Ь =-,а =1,а =!,а =1,а = —,а =-1 2! 3 2 8 3 В 4 В' ! В' 4.3 ' 3.2 ' 4.! ' 3.! 3' 4,2 Подбор параметра для интегрирования Гамильтоновых систем Для однопараметрического двухстадийного метода можно распорядиться свободным параметром так, чтобы выполнялось какое-нибудь дополнительное условие. Рассмотрим для примера случай Гамильтоновой системы, для которой движения происходят на интегральной поверхности полной энергии системы. Эта поверхность, или энергетический уровень, специфицируется заданием начальных значений координат и импульсов (точки в фазовом пространстве).
Известно, что использование методов Рунге-Кутты приводит к тому, что в процессе вычислений полная энергия может меняться. Для демонстрации эффекта различного выбора свободного параметра на поддержание постоянства энергии рассмотрим простую систему второго порядка — ангармонический маятник.
Запишем гамильтониан системы: 2 2 4 уг у, у, > гепап; Нат:ау-» — + — +— 2 2 4 теперь получим систему дифференциальных уравнений: > 1:= оларр! — Нап2(у), — Напау) у 3 1у2' у! у! Напишем процедуру для реализации одного шага по формуле метода Рунге-Кут- ты со свободным параметром Ь2: > гл23гер >м ргос (Ь2, Ь, у) 1оса1 11, Ь2! М Ь2 = О О2сп ЕККОК(Ъ2 = О') сгк! Н; Ь1 ! Ьхт(у)! Ь2:= Ьхг(у + 411(2хЬ2) ); у+ (1- Ь2) И1+ Ь2хт2 Мгп Рого Аргументами процедуры являются параметр Ь2, величина шага Ь и вектор неипвестныху (переменная типа 11зс). Результатом является переменная'типа 11эт. Чтобйполучить решение нв некотором отрезке и нарисовать фазовую траекторию, подгвтговньа:процедуру вычисления я! пгагпв! 2бб Глава 10.
Примеры решения задач л[2:ш ртос (6, Ь, у, т) 1оса1 К уу, яа1; [Г О = 0 тьсп ЕККОК('0 = 0') спи Ы; уулиу' ю1:= уу; (ог (т ю тл т)о уу:ш т1с2мср(0, Ь, Уу); яа1:= яо1, уу епв до: [ю1) опт[ ртос Для рисования графиков изменения гамильтониана от времени введем функцию, готовящую список из пар «время/гамильтонианьс > тяег)еь:= у-+ [юЧ([1Ь, Нат(у„)),1= 1 .. поря(у)П Выберем начальные данные, число шагов, величину шага Ь и несколько значений параметра Ь2: у;= [О, .4); т:= 250; Ь .ш 1; раг и> [.2,.3, .5, .85) Проведем расчеты при заданных переменной рзг значениях параметра: 1от ) то поря(рог) т[о тН [[1: = тяспея( т)с2(рог., Л, у, т ) ) епд до 1 Теперь нарисуем графики изменения гамильтониана от времени для различных значений параметра Ь2 (см.
рнс. 10.4): > р1ой [тН [[ (1 .. поря(раг) ) ), О .. т Ь, 1екелп = гоар(салтсгт, раг, шила ), 1аЬеы = [""„"Нлт), Ш[салеяь = [яея() У = 1 ., поря( раг) П, со1аг = Ыас(с лхся = вохе) о оьоя Няп\ оопп =7 Рис. 10.4. Динамика численного решения при различных константах и«тоха Проведем расчеты при тех же значениях параметра Ь2,.выбрав начальные значения, отвечаюшие большей величине энергии (см. рис. 10.5). Приведем соответствуюшие команды в обычной нотации: > у;-[0.1): и:-250: Ь:-.1: раг:-[.2..3..5А.85); > Гог 1 то поря(рвг) со тН[[1:-тьегтеь(гх2(раг[1).Ь.у,ш)):а): > р)от([тН[[(1..поря(раг)П,О..ш«|т. )ерем[-еар(сопчегс,раг,ьсппо).)аЬе)ь-["","Нап"), ььтскпеьь-[ьеч(5.5-1..поря(раг) П .со)ог Ы)асв,ахея-[охео): Движение шарика в потенциальной яие 2бб 0.91 о пи О.о наш 9.499 О 49 О ам — ыГ Рис. 10.9. Динаиика численного решения при различных константах метода На полученных рисунках даны кривые изменения полной энергии для различных значений параметра Ь2.
Видно, что в зависимости от параметра метода и энергетического уровня в расчетах получается как рост, так и убывание энергии. Для выбранного энергетического уровня и заданного шага интегрирования можно подобрать значение параметра так, чтобы в расчетах происходили оснилляпии гамильтониана около его начального значения. Движение шарика в потенциальной яме Продемонстрируем возможности пакета Мар)е для исследования обыкновенных дифференциальных уравнений и визуализации результатов на примере системы уравнений, описывающих движения шарика единичной массы.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
