Кирьянов. Самоучитель MathCad 11. Кирьянов д в
 Скачать 10.75 Mb.
|
ОбыкновенныеКутты, будет рассчитывать решение дифференциального уравнения. Чем больше step, тем с лучшей точностью будет получен результат, но тем больше времени будет затрачено на его поиск. Помните, что подбором этого параметра можно заметно (в несколько раз) ускорить расчеты без существенного ухудшения Пример решения задачи Коши для ОДУ первого порядка посредством вычислительного блока приведен в листинге Листинг 11.1. Решение задачи Коши для ОДУ первого порядка i v e n у у f t ) ( t ) d t = 0 . у := O d e s o l v e , 10 ) забывайте о том, что вставлять логические операторы следует при помощи панели инструментов Boolean (Булевы операторы. При вводе с клавиатуры помните, что логическому знаку равенства соответствует сочетание клавиш иначе будет выдано сообщение об ошибке. Как можно заметить, результатом применения блока является функция у определенная на промежутке Следует воспользоваться обычными средствами Mathcad, чтобы построить ее график или получить значение функции в какой-либо точке указанного интервала, например Пользователь имеет возможность выбирать между двумя модификациями численного метода Рунге-Кутты. Для смены метода необходимо нажатием правой кнопки мыши на области функции odesolve вызвать контекстное меню и выбрать в нем один из двух пунктов Fixed (Фиксированный шаг) или Adaptive (Адаптивный. По умолчанию применяется первый из них, т. е. метод Рунге-Кутты с фиксированным шагом Подробнее о различии этих методов сказано в разд 270 Часть III. Численные методы Встроенные функции Альтернативный метод решения ОДУ перешел из прежних версий Он заключается в использовании одной из встроенных функций rkfixed, Rkadapt или Bulstoer. Этот способ несколько первому ив простоте, ив наглядности. Поэтому я советую предпочесть вычислительный блок Однако иногда приходится решать ОДУ первого порядка с помощью второго способа, например, при следующих обстоятельствах одно ОДУ решается в контексте решения более сложных задач, в которые входят системы дифференциальных уравнений (для которых вычислительный блок неприменим в этом случае может потребоваться единый стиль программирования ответ предпочтительнее получить в виде вектора, а не функции Вы привыкли к записи ОДУ в старых версиях Mathcad, у Вас много документов, созданных сих помощью и т. п. Поскольку решение вторым способом одного ОДУ мало чем отличается отрешения систем ОДУ (см. разд приведем пример его использования в задаче из листинга 11.1 практически без комментариев (см. листинги с помощью одной из трех существующих для этих целей встроенных функций rkfixed. Обратите внимание только на необходимость явного задания количества точек интегрирования ОДУ в третьей строке листинга, а также на получение результата, в отличие от вычислительного блока, не в виде функции, а в виде матрицы размерности Она состоит из двух столбцов водном находятся значения аргумента t (от to до включительно, а в другом соответствующие значения искомой функции Листинг 11.2. Решение задачи Коши для ОДУ первого порядка вторым способом у := 0.1 D , ум := 100 у Совет В листинге 11.2 приведен пример не лучшего стиля вания. Сначала переменной у присвоено значение скаляра а затем этой же переменной присвоено матричное значение (результат решения ОДУ. Старайтесь избегать такого стиля, который ухудшает читаемость программы и может приводить, в более сложных случаях, к трудно опознаваемым ошибкам. Неплохим решением было бы назвать результат по-другому, например и Глава Обыкновенные дифференциальные уравнения 271 График решения рассматриваемого уравнения показан на рис Обратите внимание, что он соответствует получению решения в матричном виде (листинг 11.2), поэтому по осям отложены соответствующие столбцы, выделенные из матрицы у оператором <>. Примечание Пример, решенный в листингах взят из области математической экологии и описывает динамику популяций с внутривидовой конкуренцией. Сначала происходит рост численности популяции, близкий к экспоненциальному, а затем выход на стационарное состояние. Рис. Решение уравнения (листинг ОДУ высшего порядка Обыкновенное дифференциальное уравнение с неизвестной функцией y ( t ) в которое входят производные этой функции вплоть до называется ОДУ ГО порядка. Если имеется такое уравнение, то для корректной постановки задачи требуется задать начальных условий на саму функцию) и ее производные от первого до порядка включительно. В Mathcad можно решать ОДУ высших порядков как с помощью вычислительного блока таки путем сведения их к системам уравнений первого порядка. Внутри вычислительного блока ОДУ должно быть линейно относительно старшей производной, т. е. фактически должно быть поставлено в стандартной форме; • начальные условия должны иметь форму или а небо- лее сложную (как, например, встречающаяся в некоторых математических приложениях форма у Часть III. Численные методы В остальном, решение ОДУ высших порядков ничем не отличается отрешения уравнений первого порядка (см. разд. 11.1), что иллюстрируется листингом. Как Вы помните, допустимо написание производной как в виде знака дифференциала (так в листинге введено само уравнение, таки с помощью штриха (так введено начальное условие для первой производной. Не забывайте пользоваться булевыми операторами при вводе уравнений и начальных условий. Полученное решение y(t) показано на рис. Листинг Решение задачи Коши для ОДУ второго порядка = у' (0) = у := ( t , 50 Решение уравнения осциллятора (листинг Примечание В листинге решено уравнение затухающего гармонического осциллятора, которое описывает, например, колебания маятника. Для модели маятника y ( t описывает изменения угла его отклонения от вертикали, y ' ( t ) — угловую скорость маятника, y " ( t ) а начальные условия, соответственно, начальное отклонение маятника у и начальную скорость у = Второй способ решения ОДУ высшего порядка связан со сведением его к эквивалентной системе ОДУ первого порядка. Покажем на том же примере из листинга как это делается. Действительно, если формально обозна- Глава 11. Обыкновенные дифференциальные уравнения 273 чить a ( t ) , то исходное уравнение запишется через функции ив виде системы двух ОДУ • + 1 • Именно эта система решается в качестве примера в разд. 11.3. Таким образом, любое ОДУ ГО порядка, линейное относительно высшей производной, можно свести к эквивалентной системе N дифференциальных уравнений Системы ОДУ первого порядка требует, чтобы система дифференциальных уравнений была представлена в стандартной форме , t ) , , . . . , t ) О , t Задание системы (1) эквивалентно следующему векторному представлению: У < t ) = F < Y ( t > ( 2 где Y и — соответствующие неизвестные векторные функции переменной t размера векторная функция того же размера и {N+1} количества переменных (N компонент вектора и, возможно, t). Именно векторное представление (2) используется для ввода системы ОДУ в среде Для того чтобы определить задачу для системы ОДУ, следует определить еще начальных условий, задающих значение каждой из функций в начальной точке интегрирования системы to. В векторной форме они могут быть записаны в виде (где в — вектор начальных условий размера составленный из Примечание Как Вы заметили, задача сформулирована для систем ОДУ первого порядка. Однако если в систему входят и уравнения высших порядков, то ее можно свести к системе большего числа уравнений первого порядка, подобно тому как это было сделано на примере уравнения осциллятора (см. разд. Обратите внимание на необходимость векторной записи как самого уравнения, таки начального условия. В случае одного ОДУ соответствующие векторы имеют только один элемента в случае системы уравнений — N. Часть III. Численные методы Встроенные функции для решения систем ОДУ В Mathcad имеются три встроенные функции, которые позволяют решать поставленную в форме задачу Коши различными численными методами метод Рунге-Кутты с фиксированным шагом t o , — метод Рунге-Кутты с переменным шагом — метод — вектор начальных значений в точке to размера — начальная точка расчета — конечная точка расчета — число шагов, на которых численный метод находит решение — векторная функция размера двух аргументов — скалярного t и векторного у. При этому искомая векторная функция аргумента t того же размера Внимание! Соблюдайте регистр первой буквы рассматриваемых функций, поскольку это влияет на выбор алгоритма счета, в отличие от многих других встроенных функций Mathcad, например (см. разд. Каждая из приведенных функций выдает решение в виде матрицы размера В ее левом столбце находятся значения аргумента t, делящие интервал на равномерные шаги, а в остальных столбцах — значения искомых функций рассчитанные для этих значений аргумента. Поскольку всего точек (помимо начальной то строк в матрице решения будет всего В подавляющем большинстве случаев достаточно использовать первую функцию rkfixed, как это показано и листинге на примере решения системы ОДУ модели осциллятора с затуханием (см. разд. Листинг 11.4. Решение системы двух ОДУ ( t , у ) := -УО - 0 . 1 := 1 0 и ( у О , О , Глава Обыкновенные дифференциальные уравнения 275 Самая важная — это первая строка листинга, в которой, собственно, определяется система ОДУ. Сравните рассматриваемую систему (разд. записанную в стандартной форме, с формальной ее записью в чтобы не делать впоследствии ошибок. Во-первых, функция D, входящая в число параметров встроенных функций для решения ОДУ, должна быть функцией обязательно двух аргументов. Во-вторых, второй ее аргумент должен быть вектором того же размера, что и сама функция в. В-третьих, точно такой же размер должен быть и у вектора начальных значений уо (он определен во второй строке листинга). Не забывайте, что векторную функцию) следует определять через компоненты вектора ус помощью кнопки нижнего индекса (Subscript) с наборной панели Calculator (Калькулятор) или нажатием клавиши <[>. В третьей строке листинга определено число шагов, на которых рассчитывается решение, а его последняя строка присваивает матричной переменной и результат действия функции Решение системы ОДУ будет осуществлено на промежутке Как выглядит все решение, показано на рис. 11.3. Размер полученной матрицы будет равен Т е. Просмотреть все компоненты матрицы которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки. Как нетрудно сообразить, на этом рисунке отмечено выделением расчетное значение первого искомого вектора нам шаге Это соответствует, с математической точки зрения, найденному значению .07. Для вывода элементов решения в последней точке интервала используйте выражения типа .523хю 3 Внимание! -I у 2 а •4 7 га и 3 35 5 55 6 65 7 75 01 0 056 -0 -0 -0 В 057 -0 021 В 051 ЮМ 0 071 0 056 -0 047 -0 -0 -0 082 -0 -0 013 0 0 078 0 075 0 054 0 021 А • 11.3. Матрица решений системы уравнений (листинг Обратите внимание на некоторое разночтение в обозначении индексов вектора начальных условий и матрицы решения. В ее первом столбце собраны значения нулевой компоненты искомого вектора, во втором первой компоненты и т. д. Чтобы построить график решения, надо отложить соответствующие компоненты матрицы решения по координатным осям значения аргумента вдоль оси ха и — вдоль оси у (рис Как известно, решения обыкновенных дифференциальных уравнений часто удобнее изображать не Часть Численные методы в таком виде, а в фазовом пространстве, по каждой из осей которого откладываются значения каждой из найденных функций. При этом аргумент входит в них лишь параметрически. В рассматриваемом случае двух ОДУ такой график — фазовый портрет системы — является кривой на фазовой плоскости и поэтому особенно нагляден. Он изображен на 11.5 (слева, и можно заметить, что для его построения потребовалось лишь поменять метки осей на и соответственно. Примечание Фазовый портрет типа, изображенного на рис имеет одну стационарную точку (аттрактор на которую "накручивается" решение. В теории динамических систем аттрактор такого типа называется фокусом График решения системы ОДУ (листинг 0.1 11.5. Фазовый портрет решения системы ОДУ при (слева) и (справа) (листинг В общем случае, если система состоит из ОДУ, то фазовое пространство является мерным. При N>3 наглядность теряется, и для визуализации фазового портрета приходится строить его различные проекции Глава Обыкновенные дифференциальные уравнения На том же рис справа, показан для сравнения результат расчета фазового портрета с большим числом шагов. Видно, что в этом случае обеспечивается лучшая точность ив результате, решение получается более гладким. Конечно, при этом увеличивается и время расчетов. Внимание! При решении систем ОДУ многие проблемы могут быть устранены при помощи простой попытки увеличить число шагов численного метода. В частности, сделайте так при неожиданном возникновении ошибки "Found a number with a magnitude greater than (Найдено число, превышающее значение. Данная ошибка может означать не то, что решение в действительности расходится, а просто недостаток шагов для корректной работы численного ал- горитма. В заключение следует сказать несколько слов об особенностях различных численных методов. Все они основаны на аппроксимации дифференциальных уравнений разностными аналогами. В зависимости от конкретной формы аппроксимации, получаются алгоритмы различной точности и быстродействия. В Mathcad использован наиболее популярный алгоритм Рунге- Кутты четвертого порядка, описанный в большинстве книг по методам вычислений. Он обеспечивает малую погрешность для широкого класса систем ОДУ за исключением жестких систем. Поэтому в большинстве случаев стоит применять функцию rkfixed. Если по различным причинам время расчетов становится критичным или точность неудовлетворительна, стоит попробовать вместо другие функции, благо сделать это очень просто, благодаря одинаковому набору параметров. Для этого нужно только поменять имя функции в программе. Функция Rkadapt может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо, либо существуют участки медленных и быстрых его изменений. Метод с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений — частыми. В результате, для достижения одинаковой точности требуется меньшее число шагов, чем для rkfixed. Метод Булирша-Штера часто оказывается более эффективным для поиска гладких решений Решение систем ОДУ в одной заданной точке Зачастую при решении дифференциальных уравнений требуется определить значения искомых функций не на всем интервале а только водной его последней точке. Например, весьма распространены задачи поиска аттракторов динамических систем. Известно, что для широкого класса ОДУ одна и та же система при разных (или даже любых, как рассмотренный выше пример осциллятора с затуханием) начальных условиях при прихо- 278 Часть III. Численные методы дит в одну и туже точку (аттрактор. Поэтому часто нужно определить именно эту точку. Такая задача требует меньше ресурсов компьютера, чем решение системы ОДУ на всем интервале, поэтому в имеются модификации встроенных функций и Bulstoer. Они имеют несколько другой набор параметров и работают быстрее своих аналогов rkadapt to, s) — метод с переменным шагом метод — вектор начальных значений в точке to; • — начальная и конечная точки расчета — погрешность вычисления (чем она меньше, тем с лучшей точностью будет найдено решение рекомендуется выбирать значения погрешности в районе — векторная функция, задающая систему ОДУ; • к — максимальное число шагов, на которых численный метод будет находить решение — минимально допустимая величина шага. Как легко заметить, вместо числа шагов на интервале интегрирования ОДУ, в этих функциях необходимо задать точность расчета численным методом значения функций в последней точке. В этом смысле параметр похож на константу которая влияет на большинство встроенных алгоритмов Mathcad. Количество шагов и их расположение определяется численным методом автоматически, чтобы обеспечить эту точность. Два последних параметра нужны для того, чтобы пользователь мог искусственно повлиять на разбиение интервала на шаги. Параметр к служит для того, чтобы шагов не было чрезмерно много, причем, нельзя сделать к>юоо. Параметр для того чтобы ни один шаг не был слишком малым для появления больших погрешностей при разностной аппроксимации дифференциальных уравнений внутри алгоритма. Эти параметры следует задавать явно, исходя из свойств конкретной системы ОДУ. Как правило, проведя ряд тестовых расчетов, можно подобрать их оптимальный набор для каждого конкретного случая. Пример использования функции b u l s t o e r для того же примера (приведен в листинге В его первых двух строках, как обычно, определяется система уравнений и начальные условия в следующей строке матрице присваивается решение, полученное с помощью Структура этой матрицы в точности такая же, как ив случае решения системы ОДУ посредством уже рассмотренных нами ранее встроенных Глава 11. Обыкновенные дифференциальные уравнения 279 функций (см. разд. 11.3.1). Однако в данном случае нам интересна только последняя точка интервала. Поскольку сделанное численным методом количество шагов, те. размер матрицы и, заранее неизвестно, то его необходимо предварительно определить. Это сделано в четвертой строке листинга, присваивающей это число переменной в этой же строке оно выведено на экран. В предпоследней строке листинга осуществлен вывод решения системы ОДУ на конце интервала, те. в точке t 50 в виде вектора. В последней строке для примера еще раз выводится искомое значение первой функции из системы ОДУ (сравните его с соответствующим местом вектора из предыдущей строки). Листинг 11.5. Решение системы двух ОДУ , у) У - |