Методы Рунге-Кутты третьего и четвертого порядков

Рассмотрим две различные схемы Рунге-Кутты, предназначенные для численного решения обыкновенных дифференциальных уравнений первого порядка и имеющие третий порядок аппроксимации:

И две схемы Рунге-Кутты, имеющие четвертый порядок аппроксимации:

Пример. Решить методом Рунге-Кутты четвертого порядка уравнение d y /d x = – y, y (0) = 1.

В соответствии с приведенными выше соотношениями определяем коэффициенты:

Построим последовательность значений искомой функции:

…

Результаты получаемого численного решения для значения аргумента x = 10 при различных шагах интегрирования приведены в табл. 15.2. Три верные значащие цифры получены для шага h = 0.25.

Таблица 15.2. Результаты численного решения yn методом Рунге-Кутты четвертого порядка дифференциального уравнения y ' = – y с начальным условием y (0) = 1

Величина шага h 0.5 0.25 0.1 0.01 0.001 0.0001 Число шагов n         10 000 100 000 yn · 104 0.457608 0.454181 0.454003 0.453999 0.453999 0.453999

Сравнение таблиц 15.1 и 15.2 с решениями одной и той же задачи позволяет сделать вывод, что более высокая степень аппроксимации дифференциального уравнения разностным аналогом позволяет получать более точное решение при более крупном шаге и, следовательно, меньшем числе шагов, то есть приводит к снижению требуемых ресурсов ЭВМ.

На сегодняшний день для грубого расчета вычисления производятся методом Эйлера, для точного расчета — методом Рунге-Кутты.

Методы прогноза и коррекции
(итерационные методы)

Изученные нами ранее методы обладали одной важной особенностью — каждому методу соответствует обычно определенный класс точности, который мы обозначали как Oⁱ. Например, метод Эйлера обладал первым классом точности O ¹. Это означало, что с уменьшением шага в 10 раз (на порядок) точность результата повышается тоже в 10 раз (на один порядок). Метод Рунге-Кутты обладает 4 порядком точности — O ⁴, при уменьшении шага в 10 раз, результат улучшается в 10 000 раз. Поскольку этот метод по сравнению с методом Эйлера использует всего в 4 раза больше вычислений, то использование его более выгодно. На сегодняшний день известны методы до 8 порядка точности (например, метод Prince Dortmund), хотя одновременно стоит иметь в виду, что написание алгоритмов для них — задача достаточно трудная. Достоинством всех этих алгоритмов является то, что объем вычислений для них заранее известен.

Если требуется достичь ЛЮБОЙ точности на шаге, то следует использовать методы прогноза и коррекции. Этот подход состоит в том, что расчет траектории, задаваемой уравнением, на каждом шаге происходит многократно. А именно, сначала происходит расчет приближенного значения функции на конце шага какой-либо простой формулой (например, методом Эйлера), далее в этой точке вычисляется производная, и расчет происходит снова из начальной точки на шаге, но с уточненным значением производной. Последняя операция — уточнения производной и значения функции на конце шага — происходит МНОГОКРАТНО НА КАЖДОМ ШАГЕ, то есть до тех пор, пока вычисленные значения (функции и производной в конце шага) не перестанут меняться или будут меняться уже незначительно, меньше чем задаваемая заранее величина ε. Только тогда можно сказать, что точность ε достигнута.

Итак, за счет итерационной процедуры на каждом отдельном шаге можно достичь любой, наперед заданной точности ε. За такое достоинство метода приходится платить: к сожалению, невозможно сказать заранее, сколько итераций потребуется для достижения на шаге заданной точности ε. Поэтому такие методы нельзя, например, использовать в системах реального времени.

Рассмотрим для примера два метода из этого класса. Как и ранее задача состоит в нахождении функции y (t) из дифференциального уравнения dy / dt = f (y, x, t) или множества функций из системы таких уравнений.

Метод Эйлера с итерациями

1) Предсказывающая формула вычисляет (прогнозирует) значение функции на правом конце шага: y_{k + 1} = y_k + f_k · Δ t.

2) Расчитывается производная в точке k + 1 подстановкой y в исходное уравнение в k + 1 точке: f_{k + 1} = f (t + Δ t, y_{k + 1}).

3) Уточняющая формула, используя старое значение производной (с шага 1) и уточненное с шага 2, дает уточненное значение y_{k + 1}: y_{k + 1} = y_k + (f_k + f_{k + 1}) · Δ t /2. Здесь же производится подсчет итераций счетчиком i: i:= i + 1.

4) Проверка точности: | y_{k + 1} ^{i -я итерация} – y_{k + 1}^{(i + 1)-я итерация}| ≤ ε. Если условие выполнено и точность ε достигнута, то переходим на следующий шаг 5), иначе осуществляется переход на шаг 2) и процесс уточнения повторяется с новыми значениями y и f, причем их старое значение берется с предыдущей итерации.

5) Подготовка к новому шагу: изменение счетчика времени t на величину Δ t и изменение номера шага k:
t:= t + Δ t
k:= k + 1.

6) Проверка окончания расчета: t ≤ T. Если условие выполняется, расчет продолжается для следующей точки, переход на 1), иначе — конец.

Метод Милна

1) По предсказывающей формуле вычисляется грубое значение y на правом конце интервала: y_{k + 1}: y_{k + 1} = y_{k – 3} + 4/3 · (2 · f_k – f_{k – 1} + 2 · f_{k – 2}) · Δ t.

2) Рассчитывается производная в k + 1 точке: f_{k + 1} = f (t + Δ t, y_{k + 1}).

3) Снова рассчитывается y_{k + 1} по уточненной формуле, используя уже новое значение производной в точке k + 1: y_{k + 1} = y_{k – 1} + 1/3 · (f_{k + 1} + 4 · f_k + f_{k – 1}) · Δ t.

4) Рассчитывается производная в k + 1 точке с учетом вновь вычисленного более точного значения y_{k + 1}: f_{k + 1} = f (t + Δ t, y_{k + 1}). Здесь же производится подсчет итераций счетчиком i: i:= i + 1.

5) Проверка точности: | y_{k + 1} ^{i -я итерация} – y_{k + 1}^{(i + 1)-я итерация}| ≤ ε. Если условие выполнено, и точность ε достигнута, то переходим на следующий шаг 6, иначе осуществляется переход на шаг 3 и процесс уточнения повторяется с новыми значениями y и f, причем их старое значение берется с предыдущей итерации.

6) Подготовка к новому шагу: изменение счетчика времени t, изменение номера шага k:
t:= t + Δ t
k:= k + 1.

7) Проверка окончания расчета: t ≤ T. Если условие выполняется, то расчет продолжается для следующей точки, и осуществляется переход на шаг 1, иначе — конец.

Моделирование систем
с распределенными параметрами

Если объект характеризуется некоторым параметром, различным по своему значению в разных точках объекта, то можно сказать, что значения такого параметра распределены (по объекту). Если таких параметров несколько, то объект рассматривается как система с распределенными параметрами. Для осуществления расчетов систему в таком случае удобно разбить на элементарные объемы (слои). Покажем это на примере.

Пример 1. Рассмотрим процесс сушки материала (см. рис. 17.1). Допустим, что материал представляет собой кучу сырья, над которой установлен нагреватель. В процессе сушки тепло неравномерно проникает в глубь кучи, и в различных слоях материал имеет различную температуру. То есть в данном примере можно сказать, что параметр температуры распределен по объекту «куча». Опишем каждый слой своим уравнением, используя процедуру построения модели для отдельного слоя из лекции 11.

Рис. 17.1. Схема технологического процесса сушки материала, разбитого на слои

Изменение температуры в первом слое (T_н — температура нагревателя, T_ср — температура среды):

Изменение температуры в i-ом слое:

Изменение температуры в n-ом слое:

Тогда поведение системы «кучи сырья» опишется системой дифференциальных уравнений, каждое из которых опишет отдельный слой «кучи».

Алгоритм расчета такой системы показан на рис. 17.2. Особенностью алгоритма является то, что он содержит в дополнение к циклу по времени (см. алгоритм на рис. 10.5) вложенный цикл по номеру слоя. Действительно, на каждом такте необходимо отдельно просчитать изменения в каждом из слоев «кучи». Распределенность объекта (введение дополнительной координаты) имитируется дополнительным циклическим блоком.

Рис. 17.2. Блок-схема моделирования распределенной системы (на примере процесса сушки)

Опять следует обратить внимание на то, что цикл расчета производных и цикл расчета нового состояния системы отделены друг от друга во избежание появления «эффекта гонок».

Итак, инструментом моделирования распределенных систем являются вложенные циклы, по крайней мере, двойные — внутри цикла «по времени» содержится цикл «по пространству».

Модель может выглядеть по-разному — в зависимости от платформы, на которой она реализуется. Выше (см. рис. 17.2) был разобран пример, показывающий алгоритмический способ описания расчета поведения распределенной системы, использующий язык и архитектуру фон неймановской машины (см. рис. 17.3).

Рис. 17.3. Фон неймановская архитектура вычислительной системы (с общей шиной, звезда)

Представьте, что мы используем для моделирования структуру не фон неймановской машины (см. лекцию 02 ИИ), которая состоит из МНОЖЕСТВА элементарных одинаковых вычислительных элементов, соединенных между собой в сеть (см. рис. 17.4). При такой архитектуре элементы существуют одновременно, то есть параллельно синхронно во времени и в конце каждого такта обмениваются выходными (входными) сигналами между собой, соседями по сети. Такая парадигма вычислительной платформы имитирует более адекватно окружающий нас мир.

Рис. 17.4. Не фон неймановская архитектура вычислительной системы (однородная, параллельная сеть вычислительных элементов)

Используя машину фон Неймана, архитектуру с общей шиной (звезду), мы тем самым организуем процесс вычислений алгоритмически, представляя его как последовательность шагов. Сам алгоритм выполняется на одном и том же ресурсе — единственном процессоре, который через шину (центральный узел звезды) запрашивает шаг за шагом память по поводу состояния каждого из рассчитываемых им элементов. Центральный узел (шина), таким образом, является узким местом. Процессор обрабатывает такт за тактом каждый элемент имитируемой системы, которые выстраиваются в очередь к нему на обработку.

В окружающем нас мире, который мы пытаемся имитировать, в действительности сущности автономны, живут как бы сами по себе параллельно во времени, проявляя свои свойства посредством взаимодействия с соседями по пространству.

Поставим во взаимно однозначное соответствие элементарные сущности имитируемого сложного объекта из реального мира элементарным вычислителям — каждый вычислительный элемент будет имитировать одну определенную сущность (см. рис. 17.5, а). Для этого в вычислительном элементе запишем закон функционирования именно той сущности, которую он будет имитировать (см. рис. 17.5, б).

Рис. 17.5. Иллюстрация принципа объектного отражения мира в информационных моделях

Для нашей задачи расчета теплопроводности запись закона изменения температуры для элементарного объема будет иметь следующий вид:

s:= s + (al · (sl – s) + ar · (sr – s) + au · (su – s) + ab · (sn – s)) · dt.

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

Далее свяжем вычислительные элементы связями так же, как контактируют между собой элементы реального мира.

Тогда элемент будет передавать сведения о своей температуре s левому элементу в его переменную sr, правому элементу — в его переменную sl, верхнему элементу — в его переменную sb, нижнему от него элементу — в его переменную st (см. рис. 17.5, б). at, ar, al, ab — коэффициенты теплопроводности в соответствующих направлениях. Переменная s изменяется на величину подходящих с этих направлений тепловых потоков за такт времени dt: s:= s + (...) · dt. Установив, таким образом, связи всех элементов, мы получим подобие вычислительной сети реальному распределенному объекту.

Обратите внимание, алгоритм, показанный нами на рис. 17.2, исчезает, уступая место модели. В записи модели ничто не напоминает нам о специфике среды реализации. Здесь нет алгоритмических структур, соответственно и нет необходимости в посторонних по отношению к имитируемому объекту конструкциях (if, while и т.п.). Это происходит за счет того, что мы используем более полное подобие имитируемой системы платформе, на которой производится ее моделирование. При использовании фон неймановской структуры мы вынуждены осуществлять особый дополнительный перевод с языка объекта на язык машины (процедура алгоритмизации), не имеющий отношения к реальному объекту, но имеющий отношение к среде реализации. Поэтому надо стараться уделять внимание АДЕКВАТНОСТИ платформы (архитектуры, языка описания), на которой производится моделирование, реальному устройству окружающего нас мира. Стараться установить наибольшее ПОДОБИЕ между ними еще на стадии разработки.

Итак, есть две принципиально разные вычислительные архитектуры — фон неймановская и не фон неймановская (обе предложены фон Нейманом):

• фон неймановская — последовательная машина, удобная для описания процессов, последовательностей, развертывающихся во времени, и использующая нотацию алгоритмов, задач. Это централизованная система, реализующаяся по схеме «звезда»;
• не фон неймановская — параллельная машина, удобная для описания сущностей, использующая нотацию объектов, моделей. Это децентрализованная архитектура, реализующая схему «сеть».

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

Запустите последовательно три проекта: «Дом-1», «Дом-2» и «Дом-3».

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

В модели 2 проектировщик усложнил проект, установив в него дополнительные компоненты интерфейса для удобства использования. Мы видим процесс формализации в развитии. На модель накладывается интерфейс, важным здесь является то, что он в данной нотации отделен от модели. Мы можем легко заменять отдельные элементы интерфейса независимо друг от друга.

Если мы хотим уточнить модель, считая, что температуры в разных углах комнаты могут быть разными, что в действительности именно так, то нам придется разбить комнаты на более мелкие элементы, организовать их описание и установить между ними связи. Имитация станет более точной, поэтому в модели 3 вы можете уже исследовать более тонкие эффекты — попробуйте открывать и закрывать окна и двери, менять температуры отопительных приборов и окружающей среды и посмотрите на результат (см. рис. 17.6).

Рис. 17.6. Результат работы модели распределенной системы при различных условиях и внешних воздействиях

Примечание. Справедливости ради, следует отметить, что, несмотря на то, что вы видите элементы вычислительной сети, имитируете вы в конечном итоге на последовательной машине (машине процессов, алгоритмической машине) за неимением в настоящий момент у вас не фон неймановской машины. Компенсирует (маскирует) сложности алгоритмического способа мышления среда моделирования Stratum-2000, которая имитирует работу параллельной объектно-ориентированной машины на фон неймановской структуре. Но при наличии у вас не фон неймановской машины вы бы без каких-либо изменений перенесли бы в автоматическом режиме проект на сеть вычислителей. Такой перенос принято называть отражением объекта в имитационную среду.

То есть в проекте 3 мы фактически видим параллельную объектно-ориентированную систему. Структура не фон неймановской машины поддерживает объектно-ориентированный способ моделирования (см. лекцию 32).

Моделирование систем
с распределенными параметрами
при перемещении материальных масс

Еще раз усложним задачу из предыдущей лекции.

Рассмотрим процессы нагрева, сушки и одновременно перемещения неоднородной массы. Теперь (см. рис. 18.1) будем массу не только сушить, но и в процессе сушки передвигать ее по конвейеру к потребителю. Пусть длина конвейера — L, а линейная скорость его перемещения — V. Для реалистичности усовершенствуем модель, учитывая два параметра «кучи» — температуру T и влажность W каждого ее участка. Исходное состояние сырья, находящегося в бункере 1, обозначим как T ₀ и W ₀.

Рис. 18.1. Схема процесса сушки и перемещения сырья

Для отличия различных фрагментов массы не только по слоям, но и на различных участках конвейера введем деление массы на зоны. Пусть всего будет n слоев и k зон.

Время h, в течение которого лента передвигается на одну зону, составляет:

Идея моделирования состоит в том, что мы будем осуществлять при временно неподвижном конвейере нагрев и сушку (расчет температуры и влажности) по зонам и слоям, а затем перемещать сырье из зоны в зону одномоментно, имитируя перемещение конвейера.

То есть для компьютерной имитации на цифровых машинах полезно разделить процесс сушки и процесс перемещения. Так на одном такте времени Δ t имеет смысл осуществить несколько МИКРОТАКТОВ сушки a. Именно столько микротактов лента как бы стоит на месте. После этого происходит резкий сдвиг конвейера, и сырье из одной зоны попадает в следующую зону рывком. После этого снова осуществляется процесс сушки в новом положении конвейера.

Связь величин Δ t, h и a задается формулой: Δ t = h · a или так:

Далее: T _ср. — температура среды; T _Н — температура нагревателя (T _{Н j} = g · U_j); T _тр. — температура транспортерной ленты.

…

…

…

…

Тогда поведение системы «кучи сырья» опишется системой дифференциальных уравнений, каждое из которых опишет отдельный слой и зону «кучи» как по температуре, так и по влажности.

Алгоритм расчета такой системы показан на рис. 18.2. Особенностью алгоритма является то, что он содержит, в дополнение к циклу по времени и по слою (см. алгоритм на рис. 17.2), еще один вложенный цикл по номеру зоны. Действительно, на каждом такте необходимо отдельно просчитать изменения не только в каждом из слоев «кучи», но в каждой зоне.

Рис. 18.2. Блок-схема моделирования распределенного динамического процесса с перемещением масс

Итак, инструментом моделирования систем, распределенных по нескольким координатам, являются многократно вложенные циклы — внутри цикла «по времени» содержится цикл «по одной оси пространства», внутри этого цикла — еще цикл «по другой оси пространства» и так далее. Процессы, если по возможности, разделены — цикл сушки отделен от цикла движения масс, циклы процессов стоят последовательно друг за другом.

⇐ Предыдущая6789101112131415Следующая ⇒

Поделиться:

Воспользуйтесь поиском по сайту: