Математические задачи в пакете MathCAD 12
         

Глава 9.3.4. О численных методах



В завершение раздела сделаем несколько важных замечаний относительно выбора численного алгоритма решения ОДУ и задания его параметров. Они не претендуют на общность, но, надеемся, будут весьма полезны читателю, особенно в случае возникновения проблем.

Рекомендации по выбору численного алгоритма

Все численные методы решения ОДУ основаны на аппроксимации дифференциальных уравнений разностными аналогами. В зависимости от конкретной формы аппроксимации получаются алгоритмы различной точности и быстродействия. В Mathcad использован наиболее популярный алгоритм Рунге—Кутты четвертого порядка, описанный в большинстве книг по методам вычислений. Он обеспечивает малую погрешность для широкого класса систем ОДУ за исключением жестких систем. Поэтому в большинстве случаев стоит применять функцию rkfixed. Если по различным причинам время расчетов становится критичным или точность неудовлетворительна, стоит попробовать вместо rkfixed другие функции, тем более, что сделать это очень просто благодаря одинаковому набору параметров. Для этого нужно только поменять имя функции в программе.

Функция Rkadapt может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо либо существуют участки медленных и быстрых его изменений. Метод Рунге—Кутты с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений — частыми. В результате для достижения одинаковой точности требуется меньшее число шагов, чем для


rkfixed. Метод Булирша— Штера Buistoer часто оказывается более эффективным для поиска гладких решений.

Имея в виду сделанные замечания, приведем краткую сводку алгоритмов решения задач Коши для ОДУ, отмечая, какие из встроенных функций следует использовать в конкретных случаях.

  •  Для решения единственного уравнения (любого порядка) используйте вычислительный блок Given/Odesolve.
  •  Для стандартных нежестких систем используйте алгоритм Булирша— Штера Buistoer.
  •  Для систем с участками быстро и медленно меняющихся решений используйте адаптивный алгоритм Рунге—Кутты Rkadapt.
  •  В учебных целях и для решения несложных задач можно использовать алгоритм Рунге—Кутты с фиксированным шагом Rkfixed.
  •  Для получения решения в одной конечной точке интервала используйте (в зависимости от перечисленных классов задач) одну из встроенных функций с именем, начинающимся со строчной буквы.
  •  Для жестких систем (см. разд. 9.4) используйте функции Radau stiffb или stiffr, причем если вы затрудняетесь с вводом якобиана, применяйте первую из них).


Мы не будем в этой книге приводить описание численных алгоритмов, поскольку им посвящено огромное количество книг, и читатель, несомненно, найдет соответствующую литературу с подробным объяснением их работы. Тем не менее, разберем несколько важных моментов, касающихся выбора алгоритма и его параметров.

Число шагов

Все встроенные функции требуют задания для численного метода количества шагов, разбивающих интервал интегрирования ОДУ. Это очень важный параметр, непосредственно влияющий на погрешность результата и скорость расчетов. При решении систем ОДУ многие проблемы могут быть устранены при помощи простой попытки увеличить число шагов численного метода. В частности, сделайте так при неожиданном возникновении ошибки "Found a number with a magnitude greater than 10^307" (Найдено число, превышающее значение 10307). Данная ошибка не обязательно говорит о том, что решение в действительности расходится, а возникает из-за недостатка шагов для корректной работы численного алгоритма.

В качестве иллюстрации приведем результаты простых расчетов (листинг 9.6), в которых определяется пользовательская функция и(м>, представляющая решение ОДУ в зависимости от числа шагов м. На рис. 9.8 показано несколько графиков решения системы ОДУ затухающего линейного осциллятора для нескольких значений м.



Рис. 9.8. Решение системы ОДУ в зависимости от числа узлов М (продолжение листинга 9.6)



Рис. 9.9. Фазовый портрет решения системы ОДУ при М=40 и М=200 (продолжение листинга 9.6)


Как можно заметить, несмотря на, в целом, правильное решение ОДУ в узлах (оно показано кружками и квадратами), профиль решения при малых м получается неверным. Еще лучше это видно на рис. 9.9, представляющем фазовый портрет системы с разным числом шагов. Видно, что для м=200 решение получается более гладким (конечно, при этом увеличивается и время расчетов).

Листинг 9.6. Решение системы ОДУ как функции числа шагов

Погрешность алгоритмов решения ОДУ в точке

Обратимся теперь к функциям buistoer и rkadapt (тем, что пишутся со строчной буквы), которые предназначены для нахождения решения в определенной точке интервала интегрирования системы ОДУ. В первую очередь, подчеркнем еще раз, что их нельзя применять для получения решения в промежуточных точках интервала, несмотря на их вывод в матрице-результате. На рис. 9.10 показан фазовый портрет системы ОДУ динамической модели осциллятора, полученный при помощи функции buistoer в  листинге 9.5 (см. разд. 9.3.3) и с помощью rkadapt (при соответствующей замене третьей строки листинга 9.5). Видно, что, несмотря на высокую точность, равную 10-5, и верный результат на конце интервала, график, полученный при помощи функции buistoer (рис. 9.10, а), мало напоминает правильный фазовый портрет (см. рис. 9.9), начиная быть приемлемым только при предельно допустимом для обсуждаемых функций значении асс=10-16. В данном случае конкретной системы ОДУ функция rkadapt дает качественно верное решение (рис. 9.10, б), однако его необходимая точность обеспечивается только в последней точке временного интервала.

В заключение остановимся на влиянии выбора параметра асе на расчеты. Как правило, разные численные методы работают несколько по-разному. Алгоритм Рунге—Кутты дает результат тем ближе к истинному, чем меньше выбирается параметр асе, а Булирша—Штера часто демонстрирует менее естественную зависимость у(е): даже при относительно больших е реальная погрешность остается хорошей (намного лучше метода Рунге—Кутты).



Рис. 9.10. Фазовый портрет, полученный bulstoer (a) и rkadapt (б) (продолжение листинга 9.5)


Поэтому для экономии времени расчетов в функции buistoer можно выбирать и большие асе.

Чтобы обеспечить заданную точность, данные алгоритмы могут изменять как количество шагов, разбивающих интервал (t0,t1), так и их расположение вдоль интервала. Чтобы выяснить, на сколько шагов разбивался интервал при расчетах, воспользуемся простой программой, представленной в листинге 9.7. В ней определены пользовательские функции R(е) и B(е), вычисляющие для конкретной задачи (методов Рунге—Кутты и Булирша—Штера соответственно) зависимость числа шагов (т. е. числа строк получающейся матрицы) от параметра е=асс. Графики функций R(e) и B(е) показаны на рис. 9.11. Как видно, методу Рунге—Кутты при увеличении требуемой точности требуется все возрастающее количество шагов (и, соответственно, времени на расчеты), в противоположность методу Булирша—Штера, демонстрирующему большую экономичность.

ПРИМЕЧАНИЕ

Максимальное количество шагов алгоритма ограничивается еще одним параметром k встроенных функций buistoer и rkadapt (с/и. рази 9.3.3).


Листинг 9.7. Число шагов адаптивного численного алгоритма в зависимости от погрешности acc


Таким образом, проводя тестовые расчеты для различных задач и подбирая наилучший набор параметров, можно существенно сэкономить ресурсы компьютера. Конечно, проводить подобный анализ стоит в случаях, когда время расчетов для вас критично.



Рис. 9.11. Зависимость числа шагов от параметра асе численных методов (продолжение листинга 9.7)


Содержание раздела