Говорухин В., Цибулин Б. Компьютер в математическом исследовании (1185927), страница 75
Текст из файла (страница 75)
Численный анализ в МАТ[АЗ Теперь запустим команду и нарисуем исходные и отфильтрованные данные: » а-т1)(ег(о,а.х); р)ос((,х,с,г,'о') Из-за большого шага кривая исходных данных (сплошная линия) не похожа на синусоиду, а фильтр срезает амплитуды, и получаются две синусоиды большого периода, см. рис. 15.11. -1 0 20 40 80 80 100 Рис.
15.11. действия фильтра (аеианда штет) Интегрирование дифференциальных уравнений Для решения задачи с начальными данными (задача Коши) имеется ряд функций, позволяющих интегрировать системы неавтономных дифференциальных уравнений первого порядка у' Р(ф,у) с независимой переменной 1 (обычно зто время) и вектором неизвестных функций у. В стандартном наборе МАТ[.АВ для интегрирования нежестких систем могут быть использованы функции осе23 и ос)е45, а для жестких уравнений имеются функции ос)е113, ос[е155 и обе235.
Описание используемых методов численного интегрирования задачи Коши можно найти в книге [731. Обращение ко всем функциям-интеграторам практически однотипно, поэтому для изложения будем использовать символическое имя осе зо) чег, а для примеров будем применять команду осе45 (метод Рунге — кутта порядка 4/5). итак, интегрирование уравнения или системы уравнений с правой частью, вычисляемой при помощи функции Р, производится по команде [Т.Ч)-осе во)чег("Г".Т5РАМ.ЧО.ОРТ10М5.рдйя45) Здесь Т5РАМ-[ТО ТЕ1МАЦ (или диапазон значений) — массив, определяющий интервал интегрирования (значения переменной 1, для которых будет вычислено решение), ТΠ— вектор-столбец начальных условий, а дополнительные параметры ОРТ10М5 и РАМАМ5, используемые соответственно для задания параметров и передачи параметров интегрируемой системы, могут отсутствовать Функция Р(г,у) должна быть подготовлена заранее.
Ее назначение — вычислять правую часть системы для вектора у при значении независимой переменной 1 и возвращать результат в виде вектора-столбца, Результаты решения задачи Коши размещаются в матрице У, причем каждая строка есть вектор решения, отвечающий значению независимой переменной из векто- Интегрирование дифференциальных уравнений 403 ра-столбца Т. Для того чтобы получить решения в заданные моменты времени т0, С1,..., СТ) па1, следует определить аргумент Т5РАИ-[ТО Г! ... Ст)па11: Через параметр ОРТ1ОИ5 могут быть определены значения абсолютной и относительной погрешности, максимальный шаг интегрирования и другие величины, Если аргумент ОРТ1 ОИ5 отсутствует, то интегрирование ведется при значениях, заданных по умолчанию.
Например, значение относительной погрешности Ие1 То) равно 1е- 3, и все компоненты вектора абсолютной погрешности (погрешность для каждой компоненты) АЬ5ТО) одинаковы — 1е-б, Для изменения этих и других параметров применяется функция осезег, а возможные параметры и ряд переключателей перечислены в следующих далее таблицах. Заголовок функции правой части системы дифференциальных уравнений должен быть оформлен следующим образом: Гипс!(оп сит-Е(т,у,т1аО,РАРАИ5) Здесь Т) ад — строковая переменная для отработки событий.
Заметим, что в перечне параметров функции осе45 (илн другого интегратора) должен присутствовать параметр ОРТ!ОИ5, хотя бы в виде пустого массива [], если никаких изменений параметров интегрирования не требуется. Массивы Т5РЯИ, УО и ОРТ!ОИ5 можно задавать в файле, вычисляющем правую часть Р(г.у). Тогда Т5РАИ илн УО прн обращении к о((е45 можно опустить; в этом случае интегратор перед началом работы обратится к функции Е, чтобы определить недостающие параметры: П5РАи.УО.ОРт!ОИ5] Р( Ц . П, чп1!'), «Пустые» параметры можно опустить, и тогда обращение к функции-интегратору принимает кратчайшую форму: осе ьо1чег("Р") Для создания и последующего изменения параметров интегратора имеется функ- ЦИЯ ОсезЕ!. Результатом выполнения команды будет структура ОРТ10И5: ОРТ!ОН5-опеаем "НАИЕ)", ЧАЕ!.
' НАНЕ2' ЧАь2.... ) Здесь параметрам с зарезервированными именами ИАМЕ1, ИАИЕ2,... присваиваются значения ЧАЕ!, ЧЯЕ2, .... При этом не обязательно указывать все параметры, значения не упомянутых будут присвоены по умолчанию. Более того, достаточно указать несколько первых букв, идентифицирующих имя изменяемого параметра, а регистр при этом не важен. Чтобы изменить параметры интегратора на основе существующей структуры ОЕО, следует ввести: ОРТ!ОИ5-осеаег(О(О, 'ИАИЕ!',ЧАЕОЕ1,...) Для объединения назначений старой (ОЕО) и новой (ИЕИ) структур параметров используется команда, по которой новые значения замещают значения параметров старой структуры: ОРТ!ОИ5 - опеле((Я.О.
КЕН) 404 Глава 15. Численный анализ в МАТьАВ Запуск функции ог)еве1 без параметров выводит все параметры с.присвоенными им значениями. Перечислим параметры, указав в скобках тип и значение, присвоенное по умолча- нию. Таблица 15,9. Параметры команд интегрирования дифференциальных уравнений Назначение Имя Величина относительной погрешности (скаляр, 1е-3) Величина абсолютной погрешности (скалвр или вектор, 1е-5). Когда задана скалярная погрешность, то все компоненты вектора решенив будут оцениваться одинаково, а при задании вектора — каждая компонента будет оцениваться по соответствующей компоненте вектора погрешности Козффициент уточнения (целое число, равное 4 для оие45 и 1 для остальных интеграторов) увеличивает число выводимых точек; не действует, если 1епрсл(Т5РАН)>2 Имя функции, к которой интегратор обращается после успешного шага; при запуске без параметров зто имя "опер! ос" Вектор индексов компонент вектора решения, передающихся в Оисрисгсп, по умолчанию передаются все компоненты Верхняя граница шага интегрирования (положительное число, по умолчанию задается десятая часть интервала интегрирования) Начальное значение шага интегрирования (положительное число) Максимальный порядок дпя функции о)е15з (от 1 до 5) Наличие матрицы масс (попе ( М ( М(1) / М(1,у)); по умолчанию указано 'попе', т.
е. функция с правой частью системы ОДУ при обращении Е(с,у, 'пезз') возвращает матрицу масс. Для 'М' матрица масс постоянна; для 'М(1)' матрица масс зависит от времени; для 'М(1,у)' матрица масс зависит от времени н координат йе)То! Аьзто! нет)пе Оитритгсп Ои1ри15е1 Мах5тер !Ш11 а151ер Махогиег Мдзз Иногда требуется прервать расчет прн реализации некоторого события. В этом случае свойство Ечептз должно быть обьявлено (" оп" ) и в файле правой части системы дифференциальных уравнений должна быть определена функция события, которая будет вычисляться при обращении нз интегратора Е(Т.У. ' ечеп15 ' ). В процессе расчета задачи Коши будут отслеживаться переходы через нуль для функции события и формироваться дополнительные векторы ТЕ, УЕ, 1Е. В векторе-столбце ТЕ будут храниться моменты времени, для которых произошли события, матрица УЕ составится из строк с решениями, соответствующими элементам вектора ТЕ, а вектор 1Е будет содержать информацию о том, какое событие произошло.
Функция Обер!О1 позволяет параллельно с расчетом наблюдать графики зависимости от времени всех (или указанных при помощи параметри "Оитри15е! ") компонент вектора решения; при этом должна быть задействована возможность обращения к функции вывода: ОРТ1ОН5-свезет("Оисритгсп" . 'ивер!ос ') Аналогично функции оберпа52 и ог)ерназЗ позволяют отображать в процессе расчета две или три компоненты вектора решения в виде соответствующего фазо- Интегрирование дифференциальных уравнений 405 ного портрета, а оберг1пг позволяет выводить числовую информацию о реше- нии. Имеется несколько переключателей, управляющих дополнительными возможно- стями интегратора.
Допустимы значения оп нли отт, в скобках указано значение, принятое по умолчанию. Таблица 15.10. Параметры-переключатели Имя Назначение Выведение статистики вычислительных затрат (от() Наличие матрицы Якоби (отт)1 если указано "оп", то функция с правой частью системы ОДУ при обращении е(г,у, ')асоЬ(ап') выводит матрицу Якоби СР(су Матрица Якоби не зависит от решения (отт) матрица Якоби разреженнав (отт)1 если указано "оп", то функция с правой частью системы ОДУ при обращении Г(г у, ')раггегп') выводит матрицу Якоби ОР/бу как разреженную, то есть только ненулевые элементы Правая часть системы ОДУ векгоризована (отт); если указано "оп", то при обращении Р(г,[у1 у2 ...)) выводитсв [Р(г.у1) Р(г.у2) ..) Переключатель событий (отт); если указано "оп", то функция с правой частью системы ОДУ при обращении Р(г у, 'етепгз') возвращает значения функции; см.
подробности возет()е матрица масс не зависит от времени (отт) Использование формулы дифференцирования назад в функции опе155 (отт) Оценка ошибки по норме решения (отт) 5гдгз ЗасоЬ(ап 3Сопзгапг ЗРаггегп Чесгог)зее Етепгэ МаззСопэгапг ВСЕ НогеСопгго1 Для оценки погрешности по норме погв(е)<-шах(де1То)*погв(у),АЬВТо1) нужно уста- новить переключатель "оп", по умолчанию используется более жесткий контроль точности. Узнать текущее значение параметра МАНЕ можно при помощи функции о()едег, на- пример: УА[-огвдег(ОРТ10Н5.'МАНЕ' ), при этом допускается указание первых букв имени, позволяющих идентифициро- вать параметр. Если в текущем сеансе никаких назначений не было, то возвраща- ется значение, принятое по умолчанию.
В качестве примера рассмотрим систему Лоренца и подготовим в файле)ог. в функ- цию вычисления правой части системы: типсмоп т-1ог(г,у,т)ад.з1упа,г.Ь) Г 1ог - г1дЩ Ьапг) з1г)Е ОГ [огепз ецоаг1оп х правая часть систени поренца т [зтдвв*(у(2)-У( 1)); -у(2) (г-у(3))"у( 1); -Ььу(3)+у(1)чу(2)); Зададим параметры системы и определим значения для относительной и абсолют- ной погрешности, используя команду о()езег: 406 Глава 15. Численный анализ в Р(АТ[АВ » з 10.: г 24.5;, Ь В,/3,; орс-осезеы "йе1То1 , 1е-4,'АЬзТо)'.[1е-4 1е-4 1е-5]); Вычислим решение, стартуя нз точки [ 0 1 1 ], и нарисуем график зависимости от времени второй и третьей переменных, см. рис. 15.12: » [Т.У]-оое45("1ог",[О 20),[5 7 9).орс.з,г,Ь); р101(Т.У(:.2:ЗИ 40 20 15 20 10 Рис. ЗОЛ2. Интегрирование системы Лоренца Выберем последнюю полученную точку в качестве начальной и продолжим расчет: » уО-У(епс,;), [Т.т)-(к)е45("1ог .[О 30).уО.орс,з.г.Ь): Нарисуем траекторию в фазовом пространстве системы [трехмерное пространство для системы Лоренца), см.
рис. 15.13. » р1о(3(Т(:,1) т(:.2).У(:.3)): ахтз([-10 10 -10 10 0 50)), дг(0 о(Г 40 20 0 10 10 -10 -10 Рис. 15.13. Траектория в фазовом пространстве Также имеются команды оживления фазовых картин: двумерной — соаиС и трехмерной — соя)еЬЗ. В качестве иллюстрации вычислим траекторию для системы Лоренца: » [Т.Т]-оее45("1ог".[О:.02:20).[-5 -4 19].Н.10.29.8/3): и нарисуем фазовый портрет для первой и третьей переменных, см. рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
