Как зависит ошибка усечения от шага интегрирования - Ремонт и установка крупной бытовой техники

Локальная и глобальная ошибки дискретизации при численном интегрировании дифференциальных уравнений

последовательностью точек x₀,x₁. в соответствующие моменты времени t₀,t₁. Значения точек должны удоволетворять приближенному равенству

Если специально не оговорено иное, то предполагается, что моменты времени выбираются через равные интервалы с величиной шага h>0, то есть

Аппроксимируем производную в момент времени t_k соотношением

При такой аппроксимации уравнение (1) примет вид:

Формула (2) известна как прямой метод Эйлера.

На рис.1(a) показана графическая интерпретация прямого метода Эйлера. На (k+1)-ом шаге векторное поле предполагается (локально) постоянным со значением f(x_k,t_k).

Рис.1 Иллюстрация алгоритмов (а) прямого метода Эйлера, (b) обратного метода Эйлера

Меньшее значение величины шага h в итоге дает точки аппроксимации чаще и, как демонстрирует рис.2, приводит к большей точности интегрирования, что приобретает математический смысл, поскольку (2) стремится к (1) при h->0.

Рис.2 Влияние величины шага. Уравнение dx/dt=-6x+5t -t интегрируется от x=1 прямым методом Эйлера при h=0.3 (а) и при h=0.1 (b). Точное решение показано штриховой линией.

Обратный метод Эйлера подобен прямому, но есть одно отличие в аппроксимации для производной

Такая аппроксимация дает формулу обратного метода Эйлера:

На рис.1(b) показана геометрическая интерпретация обратного метода Эйлера. На (k+1)-ом шаге векторное поле предполагается (локально) постоянным со значением f(x_k+1,t_k+1).

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

Рис.3 Та же система, что и на рис.2 проинтегрирована от x₀=1.0 с h=0.3 (a) прямым методом Эйлера, (b) обратным методом Эйлера. Точное решение показано штриховой линией.

Трапецеидальный алгоритм имеет вид:

Имеется целое семейство уравнений Рунге-Кутта второго порядка. Мы рассмотрим модифицированный алгоритм Эйлера-Коши, заданный соотношением:

Из этой формулы следует, что модифицированный алгоритм Эйлера-Коши включает два этапа. На первом этапе с помощью прямого метода Эйлера происходит перемещение на пол шага вперед к моменту времени (t_k+h/2):

На втором этапе это промежуточное значение используется для аппроксимации векторного поля с помощью итераций Эйлера прямого типа:

Как и в случае алгоритма второго порядка метод Рунге-Кутта четвертого порядка относится к явным алгоритмам. Он использует промежуточные моменты времени для для вычисления состояния в момент времени t_k+1. Следующие формулы определяют алгоритм Рунге-Кутта четвертого порядка:

Для определенных выше алгоритмов на каждом шаге требуется только одна начальная точка x_k. Такие алгоритмы называются одношаговыми. Одношаговые алгоритмы высокого порядка имеют высокую точность, но они не эффективны, если велики затраты на вычисление f. Например, алгоритм Рунге-Кутта четвертого порядка требует на шаге четыре значения f. Кроме того, на текущем шаге не выполняются оценки функций с целью их использования на последующих шагах.
В отличие от одношаговых алгоритмов, многошаговые алгоритмы повторно используют предыдущую информацию о траектории. В m-шаговом алгоритме для определения x_k+1 используют m предыдущих точек x_k, x_k-1. x_k-m+1 и значения f в этих точках. Общая формула m-шагового алгоритма имеет вид

Локальная ошибка определяется как ошибка на шаге алгоритма:

Для m-шаговых алгоритмов предполагается, что предыдущие m точек x_k-i при i=0. m-1 являются точно заданными, т.е.

Ошибка усечения — это локальная ошибка, которая получилась бы в результате выполнения алгоритма на компьютере с бесконечной точностью.
Другими словами, подразумевается, что эта локальная ошибка возникает помимо ошибки округления. Также важно помнить, что для m-шаговых алгоритмов предыдущие m точек x_k. x_k-m+1 предполагаются точно заданными.
Ошибка усечения берет свое название от алгоритмов, основанных на рядах Тейлора (например, Рунге-Кутта). Эти алгоритмы были бы точными, если бы использовались полные (бесконечные) ряды. Ошибка возникает при усечении ряда до конечного числа членов.
Ошибка усечения зависит только от алгоритма. Она не зависит от используемого компьютера и, следовательно, может быть проанализирована. Для алгоритмов Рунге-Кутта K-го порядка, при подходящих условиях, локальная ошибка усечения представляет собой

где «альфа» зависит от K, f и x_k, но не зависит от h. Для многошагового алгоритма K-го порядка локальная ошибка усечения имеет вид

Глобальнае ошибка округления — это простое накопление локальных ошибок округления. Если локальная ошибка составляет «эпсилон», то ошибка округления на единичном интервале будет

Подобно ошибке округления, локальная ошибка усечения также накапливается с каждым шагом. Для одношаговых алгоритмов K-го порядка локальная ошибка усечения составляет

Если пренебречь зависимостью a_k от x_k, то на единичном интервале времени ошибка усечения будет:

Порядок метода интегрирования.

Дата добавления: 2015-06-12 ; просмотров: 3180 ; Нарушение авторских прав

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

· ошибка дискретизации, возникающая в результате замены дифференциального уравнения (3.1) разностной аппроксимацией (3.2);

· ошибка округления, накопившаяся при выполнении арифметических операций.

При этом доминирующей является, как правило, ошибка дискретизации.

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

Введем величину , называемую глобальной ошибкой дискретизации. Отметим, что зависит от величины шага, поскольку предполагается, что значения вычисляются при определенном . Воспользуемся также стандартным обозначением для величины, стремящейся к нулю при с той же скоростью, что и . В общем случае будем говорить, что функция равна , если при величина ограничена. Проведя анализ можно показать [23], что для метода Эйлера . Это, обычно, выражают утверждением, что метод Эйлера имеет первый порядок. Практическим следствием этого факта является то, что при уменьшении приближенное решение будет все более точным, и при будет сходиться к точному решению с линейной скоростью по , т.е. при уменьшении шага вдвое ошибка дискретизации уменьшится примерно в 2 раза. Столь медленная сходимость служит препятствием использования простых методов первого порядка.

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

Рассмотрим, например, явный одношаговый метод Хьюна (или метод Рунге – Кутты второго порядка). Он определяется формулой

Сравнивая его с методом Эйлера легко заметить, что значение заменено на среднее значений функции , вычисленных в двух различных точках. Данный метод имеет ошибку дискретизации .

Наиболее знаменитым из методов Рунге – Кутты, а возможно и из всех методов численного интегрирования, является классический метод четвертого порядка, задаваемый формулой

где

;

Суть метода в том, что вектор-функция правых частей системы ОДУ определяется не только в узлах сетки, но и в промежуточных точках интервала , а величина в методе Эйлера заменена на взвешенное среднее значений , вычисленных в четырех различных точках.

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

Для многошаговых методов порядок напрямую связан с объемом информации, используемой на каждом шаге. В частности, двух и трехшаговые методы Адамса – Бишфорта, задаваемые формулами (3.10) и (3.11) имеют второй и третий порядок, соответственно.

3.2.6. Процедуры численного моделирования
с автоматическим выбором шага

Данный раздел касается не выбора того или иного метода интегрирования, а реализации самой процедуры интегрирования на ЭВМ.

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

На выбор шага, как обычно влияют два фактора – точность и устойчивость. Шаг целесообразно поддерживать таким, чтобы погрешность интегрирования не превышала допустимого значения и величина шага при этом была максимально возможной.

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

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

Первый способ реализует правило Рунге [29, 23], при котором ошибка дискретизации определяется по формуле

где – некоторая константа;

– значение фазовой переменной , полученной в точке тем же методом, что и , но только за два шага интегрирования от точки , каждый из которых в два раза меньше обычного шага ;

– порядок метода интегрирования.

Для метода Эйлера формула Рунге дает

для метода Рунге – Кутты 4-го порядка

Величина погрешности аппроксимации на каждом шаге не должна превышать допускаемой погрешности . Обычно полагают, что она составляет от 0,01 до 0,001 от текущего значения определяемой фазовой координаты .

При автоматическом выборе шага наиболее популярным является алгоритм «трех зон» [29]:

где – коэффициент изменения шага, обычно равный 2, позволяет дискретно менять шаг в процессе интегрирования. Другим вариантом является алгоритм плавного изменения шага интегрирования

где – порядок метода интегрирования.

3.2.7. Особенности численного интегрирования
технических систем

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

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

Практика показала, что наличие в технических системах объектов различной физической природы приводит к тому, что процессы в них характеризуются разнотемповостью, т.е. наличием быстрых и медленных составляющих движения. Кроме того, возможно одновременное присутствие как монотонных, так и медленно затухающих гармонических составляющих. Свойство системы содержать в своем переходном процессе существенно различные по своим спектральным характеристикам составляющие принято называть жесткостью[27].

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

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

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

Пусть исследуемая система является линейной и может быть описана матричными уравнениями состояния

Матрица называется матрицей состояния или матрицей Якоби. Собственные значения матрицы определяют устойчивость и характер переходных процессов в исследуемой системе.

Составляющие движения (обычно называемые модами), связанные с собственными значениями , лежащими в левой полуплоскости далеко от мнимой оси, соответствуют быстро протекающим и, обычно, быстро затухающим процессам в системе. Собственные значения малые по модулю, лежащие близко к мнимой оси, определяют основное движение системы.

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

Оценкой жесткости системы ОДУ обычно считается число обусловленности матрицы Якоби

где – норма матрицы .

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

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

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

Элементами матрицы Якоби являются частные производные от
нелинейной вектор функции по переменным состояния системы . Для нелинейных систем жесткость, в общем случае, не является постоянной величиной и меняется в процессе интегрирования.

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

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

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

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

Изложенные рекомендации по выбору методов интегрирования жестких систем предполагают, что исследователь хорошо знаком с особенностями объекта. Например, если речь идет о проектировании или оптимизации системы управления некоторого хорошо изученного объекта, то выбор метода интегрирования может быть проведен путем сравнения нескольких вариантов. Если же объект недостаточно исследован, то эффективными могут оказаться процедуры, обладающие элементами адаптации к особенностям объекта.

Некоторые современные моделирующие программные комплексы содержат наборы методов, расчетные формулы которых настраиваются на решаемую задачу, используя для этого оценки некоторых параметров, обычно, собственных значений якобиана. Особенно перспективными считаются явные адаптивные методы, не требующие при своей реализации вычисления матрицы Якоби и решения алгебраических уравнений [27]. Такие методы есть, в частности, среди решателей ОДУ программного комплекса «МВТУ».

3.3. Моделирование гибридных (событийно-управляемых)
технических систем

Технические системы, по определению, относятся к сложным техническим системам. Одной из особенностей этих систем является то, что поведение многих из них определяется событиями, происходящими как внутри этих систем, так и в окружающей среде. Соответственно, они обладают как непрерывной, так и дискретной динамикой, находящимися в сложном взаимодействии. Подобные системы часто называют гибридными системами [20, 21]. В отечественной литературе также используются синонимы – непрерывно-дискретные, системы с переменной структурой, реактивные, событийно-управляемые.

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

Гибридное поведение может быть связано со следующими факторами:

· совместное функционирование непрерывных и дискретных объектов. Подобное поведение является типовым для непрерывных объектов (механических, гидравлических и т.д.) управляемых цифровыми регуляторами, например, для ИСЗ. Управление может формироваться как в фиксированные моменты времени, так и асинхронно, в зависимости от изменения фазовых координат объекта;

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

· гибридное поведение, обусловленное изменением состава системы.

Все перечисленные факторы являются типичными для технических систем.

Учет дискретной динамики приводит к тому, что фазовое пространство гибридных систем разбивается на области с различным поведением, при этом фазовая траектория в зависимости от происходящих событий оказывается то в одной области фазового пространства, то в другой. Достижение фазовой траекторией границы областей является событием, приводящим к смене поведения.

События могут быть как внутренними, так и внешними. Например, при моделировании робота-манипулятора внутренние события могут быть связаны с типовыми нелинейностями механической части, а внешние – являться следствием взаимодействия со средой.

Очевидной и удобной моделью систем, управляемых событиями, является конечно-автоматная модель. Например, можно использовать конечный автомат, чтобы описать автоматическую передачу автомобиля. Передача имеет ряд состояний: парковка, нейтраль, движение, реверс и т.д. Система переходит из одного состояния в другое, когда водитель перемещает рычаг из одной позиции в другую, например, из позиции парковка в нейтральное положение.

При всей ее простоте и наглядности конечно-автоматная модель достаточно строга и формальна. Однако классическое графическое представление конечных автоматов обладает рядом недостатков. Главным недостатком является отсутствие понятия времени, что предполагает статичность состояний. Другие недостатки – отсутствие иерархии состояний, обобщения переходов, средств выражения прерываний
и продолжения нормальной работы после их обработки [11].

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

В настоящее время для моделирования систем, управляемых событиями, широко используется предложенный Д. Харелом [7] визуальный формализм – Statechart (карты состояний и переходов). Карты состояния были разработаны применительно к моделированию дискретных систем, однако они могут служить хорошей основой и для моделирования гибридных систем, т.е. с их помощью можно описать поведение системы как в дискретном, так и в непрерывном времени [8]. Основные
неграфические компоненты таких диаграмм – это событие и действие, основные графические компоненты – состояние и переход.

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

Действия –это реакции моделируемой системы на события. Подобно событиям, действия принято считать мгновенными.

Состояние –условия, в которых моделируемая система пребывает некоторое время, в течение которого она ведет себя одинаковым образом. В диаграмме переходов состояния представлены прямоугольными полями со скругленными углами.

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

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

В [25] качестве примера простой дискретной системы (частного случая гибридной системы) рассмотрена диаграмма (карта) состояний и переходов цифровых часов, представленная на рис. 3.4.

Рис. 3.4. Диаграмма состояний и переходов цифровых часов

На часах имеется две кнопки: Кнопка Режима и Кнопка Увеличения. Нажатие любой из них генерирует событие, которое может вызывать переход из одного состояния в другое. Имеются три состояния: Дисплей, Установка Часов и Установка Минут. Состояние Дисплей – начальное состояние (что обозначается стрелкой, направленной от блока перехода по умолчанию в виде черного круга). Нажатие кнопки Режимы в состоянии Дисплей вызывает появление события Режимы и переход в состояние Установка часов. В состоянии Установка Часов, событие Режимы вызывает переход к состоянию Установка Минут, тогда как событие Увеличение увеличивает текущее время (число часов), которое отображается на экране, причем это происходит без изменения состояния. Каждому состоянию часов соответствует действие, записанное ниже горизонтальной линии. Оно начинает выполняться после того, как переход в это состояние произошел.

В [21] в качестве примера рассмотрена модель, которая описывает поведение маятника в виде шарика на нити, у которого в некоторый момент времени (событие «Обрыв») рвется нить. У этой модели два состояния – «Колебания» и «Полет». Оба описываются системами дифференциальных уравнений, каждая из которых является динамической системой. Карта поведений представляет собой совокупность состояний и переходов. В любой момент времени только одно состояние является текущим.

На рис. 3.5, а показана карта состояний обрывающегося маятника, на которой кроме двух обычных состояний помещены два особых состояния – начальное и конечное.

Рис. 3.5. Карта состояний и карта поведений обрывающегося маятника

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

Рис. 3.5, а отражает качественное поведение маятника. Для получения количественной картины необходимо связать с каждым состоянием локальное действие, в частности, интегрирование систем уравнений «Модель колебаний» и «Модель полета», как это сделано на рис. 3.5, b. Карта состояний, дополненная моделями состояний, названа в [20] картой поведения гибридной системы.

Одна из доступных реализаций концепции гибридного моделирования реализована в пакете StateFlow среды MATLAB, который кратко описан в главе 5. Также следует выделить отечественный пакет Model Vision Studium, ориентированный на моделирование сложных поведений [19].

Глава 4
Автоматизированное моделирование
технических объектов

Исторический экскурс. Необходимость автоматизации процессов моделирования технических объектов возникла практически одновременно с появлением вычислительных машин. Однако, исторически, разные этапы моделирования автоматизировались в разное время. Первым этапом следует считать создание библиотек численных методов исследования систем. Сами численные методы были разработаны задолго до появления ЭВМ и предназначались, прежде всего, для решения задач небесной механики. Так как первоначально расчеты велись вручную, алгоритмы методов были хорошо отлажены и оптимизированы. К концу 70-х годов прошлого столетия были созданы специализированные коллекции численных методов практически для всех областей численного анализа.

Дальнейшие шаги на пути автоматизации моделирования были связаны с разработкой систем автоматизированного проектирования (САПР) и систем автоматизации вычислительного эксперимента – пакетов прикладных программ. Автоматизации подлежали стандартные расчеты и оформление результатов экспериментов. Как правило, эти системы создавались под определенную предметную область, прежде всего, в наукоемких отраслях (космическая, авиация и т.д.). Для создания большинства пакетов прикладных программ для численных расчетов использовался язык Фортран, хорошо приспособленный именно для этих целей. Подобные пакеты создавались годами, их модификация и развитие требовали специальных знаний в предметной области, численных методах и программировании.

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

При использовании системы автоматизированного моделирования исследователь формулирует математическую модель исследуемой системы на формальном входном языке моделирования. На ранних этапах развития САМ, когда отсутствовали возможности прямого изображения структуры модели на экране монитора, использовались методы кодирования структурной информации [1, 6]. Программа модели представляла собой последовательность вызовов процедур, каждая из которых являлась моделью структурного компонента. После связывания с модулями исполняющей системы пакета моделирования список компонентов преобразовывался в независимую исполняемую программу.

В современных системах автоматизированного моделирования, исходя из соображений удобства восприятия человеком, используются, как правило, графические средства задания исходной информации о модели. Переход к изображению структуры системы на экране монитора позволил исключить этап ручного кодирования схемы, однако потребовал разработки нового принципа организации самого процесса численных расчетов – принципа Data Flow, или принципа потока данных.

Можно отметить следующие факторы, способствующие внедрению систем автоматизированного моделирования:

· трудоемкость получения математической модели сложных технических объектов, связанная с опасностью совершить ошибку в многочисленных преобразованиях модельных выражений;

· необходимость многовариантного моделирования, при котором необходимо иметь для одного объекта несколько моделей, отличающихся по сложности;

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

4.1. Особенности современных систем
автоматизированного моделирования

Современные профессиональные САМ имеют следующие отличительные черты [31]:

· поддержка иерархического проектирования как сверху – вниз, так и снизу – вверх, за счет реализации многоуровневого моделирования и метода локальной детализации модели;

· компонентное моделирование на основе использования библиотек, содержащих большое число графических и функциональных описаний компонентов, причем эти библиотеки открыты для добавления в них новых описаний, которые может сделать сам пользователь;

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

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

· наличие в современных САПР и САМ постпроцессоров моделирования, что позволяет не только просматривать в удобной для пользователя форме результаты моделирования, но и обрабатывать эти результаты;

· наличие встроенных средств численного моделирования рабочего процесса в режиме реального времени или в режиме масштабирования модельного времени;

· реализация механизмов продвижения модельного времени, основанных как на принципе , так и на принципе ;

· интегрируемость с другими пакетами аналогичного назначения, которая обеспечивается соответствующими программами – конверторами, позволяющими импортировать и экспортировать данные из одной системы в другую;

· наличие средств, обеспечивающие формирование виртуальных аналогов измерительно-управляющей аппаратуры.

Если САМ предназначена для решения исследовательских задач, то к перечисленным качествам добавляются возможности активного вычислительного эксперимента [21]. В частности:

· визуализация результатов во время эксперимента;

· возможность интерактивного вмешательства в ход моделирования;

· возможность использования 2D и 3D анимации.

4.2. Иерархическое проектирование
и многоуровневое моделирование технических систем

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

Сложные системы имеют, как правило, иерархическую структуру. Естественные возможности человека позволяют оперативно обрабатывать не более единиц информации одновременно. В процессе проектирования исследователю удобно сосредотачиваться сначала на поведении отдельных компонентов системы, а затем на их взаимодействии [21]. При необходимости модели компонентов могут детализироваться или, наоборот, укрупняться.

Такой подход, известный как иерархическое проектирование [21, 31], является типовым при разработке сложных технических объектов и заключается в разбиении исходной задачи на подзадачи.

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

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

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

Использование принципа иерархического проектирования позволяет ограничить текущую сложность проекта на приемлемом уровне, за счет того, что в каждый момент времени разработчик имеет дело только с двумя смежными уровнями декомпозиции объекта – структурным описанием проектируемого в данный момент фрагмента и функциональным описанием внешней среды.

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

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

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

Такой возможностью обладает, в частности, пакет LabVIEW, пакет Simulink и основанные на нем пакеты SimMechanics и SimPower. Пакеты IDEF-моделирования (ERWIN, BPWIN) принципиально основаны на многоуровневом изображении объектов.

4.4. Архитектура программ
автоматизированного моделирования

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

4.4.1. Графический интерфейс программ математического
моделирования динамических систем

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

Рис. 4.1. Модульная структура
системы автоматизированного моделирования

Графический пользовательский интерфейс позволяет вводить информацию об исследуемой системе путем «рисования» на экране монитора проектируемой схемы в виде, понятном широкому кругу специалистов. Формой графического представления информации о моделируемой системе могут быть [16]:

· операторно-структурные схемы, принятые в ТАУ;

· функциональные и принципиальные схемы различных физических устройств;

· кинематические схемы механизмов;

· блок-схемы алгоритмов и другие графические модели.

Однако простым «рисованием» роль графического интерфейса не ограничивается. Задачами графического интерфейса, кроме того, могут быть [16]:

· контроль за соблюдением некоторых правил в процессе создания графического изображения на экране монитора (обычно накладываются ограничения на способы соединения компонентов и т.п.);

· преобразование информации о схеме в команды для моделирующей программы (моделятора);

· контроль за процессом моделирования, визуализация результатов моделирования.

П.5. Локальные и глобальные погрешности одношаговых методов решения ДУ

(метода Эйлера и методов Рунге-Кутта 2го, 4го порядка).

Теорема 6.1:

Если локальная погрешность метода , то глобальная .

как и при численном интегрировании, при переходе от локальной погрешности к глобальной, точность метода уменьшается на порядок. (6.8):

Методы	Локальная	Глобальная
Эйлер	const*h 2	const*h
Р.–К. 2го порядка по времени	const*h 3	const*h 2
Р.–К. 2го порядка по производной	const*h 3	const*h 2
Р.–К. 4го порядка	const*h 5	const*h 4

Как и при численном интегрировании, порядок метода – степень h в глобальной погрешности.

П.6. Многошаговые методы решения ДУ и СДУ.

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

В k-шаговом методе используются значения с k предыдущих шагов.

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

Формулы метода Милна:

(6.9)

Метод Милна – 4х шаговый (т.к. использует 4 предыдущих значения) и имеет 4-ый порядок точности. Перед применением метода Милна нам надо знать 4y, следовательно, необходимо сделать хотя бы 3 шага каким-нибудь одношаговым методом.

П.7. Оценка погрешности решения ДУ и СДУ методом двойного пересчета. Коррекция решения.

Используя такую же идею, как и в численном интегрировании, находим решение ДУ на [a,b] дважды с шагом h и с шагом h/2. Получим следующую картину:

Сравниваем попарно, если расхождение между для метода 2го порядка, для метода 4го порядка, то в качестве точного решения берём . Если же точность не достигнута, то шаг h уменьшаем вдвое и т.д., пока она не будет достигнута.

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

Как и при ЧИ, при решении ДУ и СДУ после 2го пересчёта в качестве точного решения выгодно брать не , а .

— для второго порядка

Метод двойного пересчёта применим не только лишь при ЧИ, при решении ДУ и СДУ, но и при решении других численных методов.

П.7. Краевые задачи для дифференциальных уравнений.

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

Рассмотрим простейшую краевую задачу для ДУ 2го порядка:

(6.10)

А мы умеем решать:

(6.11),

В (6.11) нам известно , поэтому для решения задачи (6.10) мы будем подбирать в (6.11), с тем, чтобы у(b) = у₁

Метод стрельб

После пристрелки и определения интервала [a,b],

где идёт смена знака, запускаем МПД или МХ.

На практике это выглядит так, как будто мы

решаем уравнение , где возвращает

решение задачи Коши (6.11) в точке b при

П.9. Что делать, если ДУ не может быть разрешено относительно старшей производной?

Так как ДУ не может быть решено относительно старшей производной, то тогда на каждом шаге решаем нелинейное уравнение относительно y ( n) (все остальные неизвестные y,y’,y”,…, y ( n-1) -к этому моменту уже известны).

Решать уравнение относительно старшей производной любым методом(Хорд, МПД, Ньютона).

Замечание:

Таким образом, если ДУ не разрешается относительно старшей производной, то у нас возникает дополнительный цикл (самый внутренний) при написании программы.

Типы и классификация ошибок численного интегрирования

Современный численный метод одновременно с решением задачи Коши должен вычислять значение или оценку ошибки. Первый такой метод был предложен Мсрсоном в 1959 году (метод Рунгс-Кутта-Мсрсона). Теоретически таких методов описано достаточно много. Однако хороших практических реализаций в виде алгоритмов или прикладных программ мало. Рассмотрим причины возникновения ошибки при решении задачи Коши. Эта ошибка складывается из следующих компонент.

1. Неустранимая погрешность исходных данных. На нее повлиять нельзя, но учесть необходимо.
2. Локальная ошибка метода. Ее величина определяется порядком р используемого метода и зависит от шага: е,

()(h p *’). Эта ошибка связана с отбрасыванием остатка ряда при разложении точного решения в ряд по шагу при условии, что вес исходные данные заданы точно и вычисления проводятся без ошибок округления.

3. Ошибка переноса, которая возникает из-за того, что ошибка, сделанная на предыдущем шаге, переносится на следующий шаг с каким- то коэффициентом.

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

Если производить расчет на малом отрезке интегрирования, то ошибкой округления можно пренебречь, и преобладающей становится ошибка метода. Но когда количество шагов велико, ошибки округления накапливаются и становятся преобладающими. Экспериментально устаО а

количестве шагов взаимодействуют, в результате возникает общая ошибка, накопленная на веем отрезке интегрирования при применении данного метода. Эта общая ошибка ?_г называется глобальной ошибкой. Известно, что при интегрировании на отрезке |0; Т|

Здесь множитель Т)(Т) зависит только от вида функции f(x) и нс зависит от шага интегрирования. Соотношение (7.9) называют правилом Рунге.

Общая ошибка ?у, порожденная всеми составляющими, сеть их сумма (рис. 7.1). Существует оптимальное значение шага /Г, обеспечивающее наименьшую суммарную ошибку. В случае, изображенном на рис. 7.1, отклонение от этого оптимального значения нс вызывает силь-

„ _ , „ , ного изменения ошибки. Возможно, однако, си-

Рнс.7.1. Суммарная ошибка

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

Полученную в результате интегрирования ошибку нужно уметь оценивать. Кроме того, полезно иметь возможность сознательно так выбирать шаг интегрирования, чтобы влиять на величину ошибки.

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

Очевидно, именно этот член и характеризует величину ошибки. Значит, следует оценить значение величины R. Будем считать, что при малом /| на протяжении шага R = const. Пусть на некотором шаге с номером г выбрана точка t е [t_r‘,t_nl]. Тогда x(t’) — точное, а х*/, — вычисленное при интегрировании с шагом /; значение искомой функции. Можно считать, что на шаге

Теперь разобьем шаг на два, длины —. и вычислим то же значение,

выполнив два шага интегрирования. Будем считать, что величина R при этом нс изменилась. Получим:

Вычитая уравнения, можно получить

И можно считать, что при интегрировании с постоянным шагом И верно соотношение

Таким образом, схема выбора величины /г на шаге для достижения заданной точности Очакова:

1. при известном х, с шагом /г вычислить т_г+/;

источники:

http://life-prog.ru/2_14935_poryadok-metoda-integrirovaniya.html

http://megalektsii.ru/s20945t2.html

http://bstudy.net/718789/ekonomika/tipy_klassifikatsiya_oshibok_chislennogo_integrirovaniya

Интегрирование функции, заданной таблично

Читайте также:

Анализ условий выполнения заданной операции, анализ опасностей и вредностей производства. Технические требования на автоматизацию операции
ВИДЫ, ФУНКЦИИ, ЛИНГВИСТИЧЕСКИЕ ОСОБЕННОСТИ
Глава 5. ИНТЕГРИРОВАНИЕ
ДАЛЬНОМЕРА С ДВОЙНЫМ ИНТЕГРИРОВАНИЕМ
Идея табличного представления множества форм. «Формации» и «категории». Двойственный подход Боулдинга к построению теории систем и его экспликация
Интеграле. Интегрирование чётных и нечётных функций.
ИНТЕГРИРОВАНИЕ НЕКОТОРЫХ ИРРАЦИОНАЛЬНОСТЕЙ

Если подынтегральная функция задана таблично в виде пар значений x(i),y(i) (узлов), то интеграл можно вычислить несколькими способами. Первый заключается в том, чтобы выразить зависимость y от x какой-либо подходящей, т. е. решить задачу аппроксимации или интерполяции табличных данных. Затем эта зависимость используется для интегрирования функции методами, описанными выше. Выбор аппроксимирующей и интерполирующей функции, а также методы расчета их параметров описаны в соответствующем разделе.

Задачу интегрирования таблично заданной функции можно решить, не прибегая к построению аппроксимирующей (интерполирующей) функции. Если табличные данные приводятся с постоянным и достаточно маленьким шагом по х, то можно применить квадратурные формулы. Пределы интегрирования могут быть любыми в пределах табличных данных и совпадать с узлами. В программе 20 выполнен расчет с помощью метода Симпсона (он оптимален для интегрирования табличных зависимостей). Единственным его недостатком является требование четности интервалов интегрирования N. А изменить количество интервалов интегрирования таблично заданной функции мы не можем. Для метода трапеций этой проблемы нет и в программе 20 приведен также расчет методом трапеций. Для иллюстрации точности интегрирования в программе 20 в качестве табличных данных взята синусоидальная зависимость с точностью до третьего знака после запятой. С такой же точностью поручены, как видно из листинга программы, и значения интегралов. Это правило выполняется всегда: чем точнее задана таблица, тем точнее можно вычислить интеграл и тем менее точный метод интегрирования можно использовать.

Если же шаг интегрирования не постоянный. А это часто бывает с экспериментальными табличными данными, то лучше воспользоваться первым способом: построить аппроксимирующую (интерполирующую) функцию. В программе 21 приведены три наиболее компактные формы интерполяции и аппроксимации табличных данных для целей интегрирования. Как видно из листинга программы, наибольшую точность дают кубическая сплайн-интерполяция и аппроксимация полиномом (в данном случае четвертой степени).

Программа 20

Программа 21

Расчет изменений термодинамических функций в ходе химической
реакции по интегральным уравнениям

В основе расчета изменений термодинамических функций: энтальпии DrH0Т, энтропии DrS0Т и энергии Гиббса DrG0Т , а также константы равновесия для химической реакции лежат соответствующие дифференциальные уравнения и их интегральные формы, представленные в таблице 7 [11].

Таблица 7

Дифференциальныe уравнения	Интегральные формы
(А)
(Б)
(В) (Г)
(Д)
Вспомогательные константы и функции

Интегральные уравнения получены на основании температурной зависимости теплоемкости индивидуального вещества:

(47)

Для химической реакции

(48)

где знак D_r обозначает изменение соответствующих коэффициентов в ходе реакции с учетом стехиометрии (иногда его называют химическим оператором), например

(49)

Все дифференциальные уравнения (А-Д) приведены к форме

(50)

Если принять T₀=298 K, то соответствующие значения y₀ можно вычислить, используя справочные данные о теплотах образования (D_rH_f⁰₂₉₈ и абсолютных энтропиях (D_rS⁰₂₉₈) веществ – участников реакции.

(51)

(52)

(53)

(54)

В программе 22 дан пример использования интегральных форм уравнений (А-Д) для расчета изменений термодинамических функций химической реакции, который может быть использован для расчетов в домашних заданиях по физической химии. В разделе «Решение дифференциальных уравнений» будут приведены программы решения дифференциальных уравнений (А-Д), а также их систем. Для расчета всех термодинамических данных, как следует из приведенных выше уравнений, достаточно рассчитать две из них, например и , так как уравнение Гиббса-Гельмгольца (53) справедливо не только при 298 К, но и при любой температуре. Поэтому для получения полной информации о термодинамики химической реакции достаточно системы из двух дифференциальных уравнений, например, А и В.

В конце программы приведен расчет изменения энергии Гиббса и логарифма константы скорости реакции при заданных температурах методом Темкина-Шварцмана. Обычно в домашних заданиях по физической химии этот расчет занимает много времени в особенности если учесть, что значения постоянных М₀, М₁, М₂ и М_-2 надо брать из справочника, в котором эти значения приведены с шагом в 100 К. Провести же расчет этих постоянных по соответствующим формулам при любых температурах в MathCad не представляет труда.

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

Программа 22

Контрольные вопросы к главе 5

1. Каким образом дифференциальное уравнение можно привести к интегральному виду?

2. Какое решение дифференциального уравнения называют «точным»?

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

4. Какие принципы лежат в основе численных методов интегрирования?

5. Как зависит ошибка усечения от шага интегрирования?

6. Какими способами можно уменьшить ошибку усечения?

7. От чего зависит ошибка округления?

8. Покажите на рис.5, где ошибка усечения в методе трапеций и в методе прямоугольников? В каком методе она меньше?

9. Какими двумя методами получить расчетные формулы метода прямоугольников и метода трапеций?

10. Поясните, как получена расчетная формула метода прямоугольников?

11. Получите расчетную формулу метода трапеций графическим методом.

12. Как получить расчетную формулу Симпсона из квадратурной формулы Котеса?

13. Как оценить точность вычислений в методах прямоугольников, трапеций и Симпсона?

14. Как задать точность вычислений определенного интеграла при использовании встроенных функций MathCad?

15. Как средствами MathCad взять интеграл и получить аналитическое решение?

16. Что обозначает термин «форматирование результата»?

17. Как отформатировать результат?

18. В чем состоит алгоритм интегрирования Ромберга?

19. Как выбрать метод интегрирования таблично заданной функции?

20. Какие численные методы можно использовать для замены табличной функции?

21. Как получены интегральные уравнения в таблице 7? Можно ли их получить, используя символьный процессор MathCad?

Расчетное многовариантное задание № 5

А. Вычислите определенный интеграл методом прямоугольников и методом трапеций с заданным количеством интервалов N; оцените точность определения интеграла каждым из этих методов;

Б. Вычислите определенный интеграл методом Симпсона с тем же количеством интервалов; сравните точность интегрирования методом Симпсона и методами трапеций и прямоугольников;

В. Вычислите определенный интеграл методом Симпсона с заданной точностью e; сравните его значение со значением, вычисленным с помощью встроенной функции MathCad; попытайтесь взять интеграл символьно.

Таблица 8

№ вар.	Интеграл	N	e
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001
			0.000001
			0.0001
			0.00001

Форма записи отчета в лабораторном журнале:

Дата:____. Занятие № __. Тема: «Вычисление определенного интеграла». Вариант ___.

а). Значение интеграла =0.23456, метод прямоуг., N= 16, достигнута точностью 1.2×10^-3

Значение интеграла =0.23387, метод трапеций , N= 16, достигнута точностью 5.2×10^-5

б). Значение интеграла =0.23388, метод Симпсона , точность 1.2×10^-6

в). Значение интеграла = 0.23389 методом Симпсона с заданной точностью 0.00001. Такое же встроенной функцией.

Имена программ: интеграл1. mcd ; интеграл2. mcd, интеграл3. mcd

Расчетное многовариантное задание № 6

А. Вычислите определенный интеграл таблично заданной функции Y=f(X) (см. данные в табл. 2). Метод выберите в зависимости от шага и количества узлов таблицы.

Б. Вычислите изменения термодинамических функций в ходе химической реакции при температурах от Т₁ до Т₂c шагом 50 К по интегральным уравнениям; постройте графики зависимости этих функций от температуры, а также зависимости и . Уравнения реакций даны в табл. 10.

Таблица 9

№ вар.	№ реакции	Т₁ , Т₂, К	№ вар.	№ реакции	Т₁ , Т₂, К
1.		700, 1000	23.		300, 500
2.		800, 1100	24.		200, 500
3.		300, 500	25.		300, 600
4.		200, 500	26.		400, 700
5.		300, 600	27.		500, 800
6.		400, 700	28.		600, 900
7.		500, 800	29.		700, 1000
8.		600, 900	30.		800, 1100
9.		700, 1000	31.		300, 500
10.		800, 1100	32.		200, 500
11.		300, 500	33.		300, 600
12.		300-500	34.		400, 700
13.		700, 1000	35.		500, 800
14.		800, 1100	36.		600, 900
15.		300, 500	37.		700, 1000
16.		200, 500	38.		800, 1100
17.		300, 600	39.		300, 500
18.		400, 700	40.		200, 500
19.		500, 800	41.		300, 600
20.		600, 900	42.		400, 700
21.		700, 1000	43.		500, 800
22.		800, 1100	44.		600, 900

Уравнения реакций:

Таблица 10

	2Н_2(г)+СО_(г) = СH₄O_(ж) _{метанол}		CO_2(г)+H_2(г)=CO_(г)+H₂O_(ж)
	4HCl_(г) + O_2(г) =2H₂O_(г)+2Cl_2(г)		CO_2(г)+4H_2(г)=CH_4(г) +2H₂O_(ж)
	NH₄Cl₍_b) = NH_3(г)+HCl_(г)		2CO₂_(г)=CO_(г)+O₂_(г)
	2N_2(г)+6H₂O_(ж) = 4NH_3(г)+3O_2(г)		CH_4(г) +CO_2(г)=2CO_(г)+2H_2(г) _метан
	4NO_(г) + 6H₂O_(Ж)=4NH_3(г)+5О_2(г)		C₂H_6(г)=C₂H_4(г) +H₂O_(ж) _{этан этилен}
	2NO₂₍_г)=2NO₍_г)+O₂_(г)		C₂H₆O_(ж) =C₂H_4(г) +H₂O_(ж) _{этанол этилен}
	N₂O₄₍_г)=2NO₂₍_г₎		C₂H₄O_(г)+H_2(г)=C₂H₆O_(ж) _{ацетальдегид этанол}
	Mg(OH)₂_кр=MgO_кр+H₂O₍_г₎		C₆H_6(ж)+3H_2(г)=C₆H_12(ж) _{бензол циклогексан}
	CaCO₃_кр= CaO_кр+CO₂₍_г₎ _{кальцит}		2H_2(г)+ CO_(г)=CH₄O_(г) _{метанол}
	Ca(OH)₂=CaO_кр+H₂O₍_г₎		2H₂O_(ж)+2Cl_2(г)=4HCl_(г)+O_2(г)
	S_ромб.+2H₂O_(ж)=SO_2(г)+2H_2(г)		S_монокл+2H₂O_(ж)=SO_2(г)+2H_2(г)
	S_ромб.+2CO_2(г)=SO_2(г)+2CO_(г)		S_монокл+2CO_2(г)=SO_2(г)+2CO_(г)
	2SO₂₍_г₎+O₂₍_г)=2SO₃₍_г)		SO_2(г)+Cl_2(г)=SO₂Cl_2(ж)
	SO_2(г)+Cl_2(г)=SO₂Cl_2(г)		C₂H₄O_(г)+H_2(г)=C₂H₆O_(ж) _{этиленоксид}_этанол
	CO_(г)+3H_2(г)=CH_4(г) +H₂O_(ж)		C₆H_6(ж)+3H_2(г)=C₆H_12(г) _бензол_{циклогексан}
	2CO_(г)+SO_2(г)=S_ромб+2CO_2(г)		2N_2(г)+6H₂O_(г) = 4NH_3(г)+3O_2(г)
	CO_(г)+Cl_2(г)=COCl_2(г)		4NO_(г) + 6H₂O_(г)=4NH_3(г)+5О_2(г)

Варианты творческих заданий

1. Модифицируйте программу 16 так, чтобы ее можно было бы использовать для расчета интеграла с заданной точностью.

2. Выполните расчеты изменений термодинамических функций в ходе химической реакции в соответствии с домашним заданием по физической химии.

Дата добавления: 2015-07-19; просмотров: 629 | Нарушение авторских прав

mybiblioteka.su — 2015-2023 год. (0.044 сек.)

Источник

будет соответствовать так называемым многошаговым формулам интегрирования.

Используя выражение (10.10), при соответствующих коэффициентах a0 ,1 и b0.1 , можно получить все три раннее рассмотренные формулы

интегрирования. Рассмотрим некоторые свойства этих формул, переписав соотношение (10.10), в виде

a x( t	0	) a x( t	0	h ) h ( b x‘ ( t	0	) b x‘ ( t	0		h ) 0 .
0				1											0					1
Разложим функции	x( t	0	h ) и	x‘ ( t	0	h ) в ряды Тейлора

a x( t		) a	[ x( t		)	h	x‘ ( t			)	h2	x» ( t			)	h3	x»’ ( t	) …]
0	0		0		0		0
0				1						1!					2!				3!

h b	x‘ ( t			) h b	[ x‘	( t		)		h	x» ( t		)	h2	x»’ ( t		) 0.
0	0			0		0
0							1									1!				2!

Перенесем слагаемые второго и более высоких порядков в правую часть и запишем

[ a a ] x( t	0	) [ a b b ] h x‘	( t	0	)
0			1				1	0		1
	a1				2	»			a1	b1	3		»’
			b1	h		x	( t0	)			h	x	( t0	) .

	2!								3!	2!

Это выражение будет удовлетворяться тождественно, при любых значениях x( t ) и ее производных в точке t t0 , в случае обращения в нуль

сомножителей в квадратных скобках. Это приводит к следующим равенствам:

1)для левой части — a0 a1 0 ; a1 b0 b1 ;

2)для правой части — a21! b1 0 ; a31! b21! 0 . Из анализа следует:

1.Выбором четырех различных коэффициентов a0 ,1 , b01 , нельзя

обратить в нуль все сомножители в квадратных скобках.

2. Наибольшее число первых сомножителей, которое может быть обращено в нуль, равно трем и это соответствует следующему выбору коэффициентов – a0 1 ; a1 1 ; b0 b1 0.5 . Подстановка их в исходное

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

значение первого множителя, отличного от нуля — c3 0.5 .

3. Если взять набор коэффициентов a0 ,1 , b0 ,1 , приводящий к формуле трапеций и принять, что искомая функция x( t ) является полиномом второй

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

Из сказанного следует, что в формуле трапеций искомая функция аппроксимирована полиномом второй степени и говорят, что метод численного интегрирования, опирающийся на метод трапеций, имеет порядок p 2 . Погрешность описания функции x( t ) является полиномом

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

В связи с этим первый, не равный нулю сомножитель в уравнении разложенным в ряд Тейлора обозначают c p 1 и называют ошибкой усечения.

Для формулы трапеций: c p 1 c3 0.5 . Определив значения коэффициентов a0 ,1 , b0 ,1 , для прямой и обратной формул Эйлера, найдем, что порядок интегрирования у этих формул равен 1 , а ошибка усечения с2

составляет, соответственно, 0.5 и 0.5 . Вспомнив, что в рассматриваемом примере, формула трапеций дала меньшую погрешность и это, как видим, связано с более высоким порядком метода и соответственно с меньшей ошибкой усечения. Более высокое значение порядка p и меньшая ошибка

усечения c p 1 предпочтительны при оценке формул численного интегрирования.

10.4. Устойчивость методов интегрирования

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

Изложим проблему устойчивости на примере простейшего дифференциального уравнения

имеющего аналитическое выражение решения, в виде x x0 exp( t ) ,

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

x( t ) Ai exp( i t ) ,

где i — корни характеристического уравнения. Таким образом, решение

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

Устойчивость прямой формулы Эйлера. Исследование устойчивости численных методов начнем с прямой формулы Эйлера

x1 x0 h x0‘ .

Используя (10.11) и, подставляя в формулу, вместо x0‘ , его значение x0 ,

получим

x1 (1 h ) x0 .

На следующем шаге интегрирования, соответственно, получим

x x h x‘ (1 h ) x (1 h )2 x .
2	1	1			1	0
Продолжив процедуру интегрирования в пределах n шагов, получим
		x	n	(1 h )n x .
				0
Предположим,	что	время	интегрирования	t и, соответственно,
n , тогда, для	того	чтобы	xn было ограниченным, для	устойчивого
дифференциального	уравнения,	когда Re 0,	необходимо	выполнение

условия

1 h 1 ,

где h — размер шага (действительное число); — корень характеристического уравнения (возможно комплексный).

Найдем области устойчивости, удовлетворяющие неравенству. Введя

обозначение

h q u j v ,

и подставив его в неравенство, получим

1 u j v 1

или

(1 u )2 v2 1 .

Это выражение описывает область устойчивости решения внутри единичного круга с центром в точке ( 1,0 ), изображенную на рисунке 10.2.

Полученный результат используется следующим образом. При Re 0 , шаг выбирается, исходя из условия q h 1 , что соответствует точке

внутри единичной окружности. В этом случае прямая формула Эйлера дает устойчивые решения, т.е. значение функции x( t ) при увеличении n будет

оставаться конечной величиной. Величина шага определяет и устойчивость, и точность численного интегрирования и из соображений устойчивости

может оказаться меньше, чем того требует точность. Так, при больших		,

шаг выбирается малым для обеспечения устойчивости.
Re 0 соответствует физически неустойчивой схеме, однако	при

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

В качестве иллюстрации, применим прямую формулу Эйлера к

простейшему дифференциальному уравнению	вида	x‘ x,	при	x 1	,
				0
соответствующему уравнению вида (10.11),	при	1 .	Результаты

численного интегрирования, при разных h , сведены в таблицу 10.4.

			75
							Таблица 10.4
Результаты интегрирования уравнения	x‘ x, при различных h , по
прямой формуле Эйлера.

	Шаг	Значения			при
	n	h=0.1		h=2		h=3
	0	1		1		1
	1	0.9		-1		-2
	2	0.81		1		4
	3	0.729		-1		-8
	4	0.65.61		1		16
	5	0.59049		-1		-32
	6	0.53144		1		64

Первая колонка решений, при h 0.1 , соответствует q h 0.1 , т.е. области внутри единичной окружности и устойчивому решению. Вторая колонка соответствует решению, при h 2 , q 2 , т.е. на границе области

устойчивости. Решение в этом случае осцилирует, но не нарастает. В третьей колонке решение, при h 3 , q 3 , находится вне области устойчивости.

Решение нарастает и осцилирует, хотя, в соответствии с аналитическим решением равным exp( t ), должно уменьшаться.

Устойчивость обратной формулы Эйлера. Продолжим исследование проблемы устойчивости для обратной формулы Эйлера

x1 x0 h x1‘ .

Используя (10.11), и подставляя в формулу, вместо x1‘ , его значение x1 ,

получим

x1 x0 h x1 ,

или

x1 x0 /(1 h ) x0 /(1 q ).

На следующем шаге интегрирования, соответственно получим x2 x1 h x‘2 x1 /(1 h ) x0 /(1 h )2 .

Продолжив процедуру интегрирования в пределах n шагов, получим xn x0 /(1 h )n x0 /(1 q )n .

Для устойчивого решения, устойчивого дифференциального уравнения, когда Re 0 , при t и, соответственно, при n , необходимо выполнение условия

1 /(1 q ) 1 .

Найдем области устойчивости, удовлетворяющие неравенству. Раскрывая обозначение q, получаем

1 (1 u )2 v2 .

Этому выражению при равенстве соответствует окружность единичного радиуса с центром в точке (1,0 ) , изображенной на рисунке 10.2. Неравенство

удовлетворяется вне окружности. Таким образом, обратная формула Эйлера устойчива для устойчивых схем, когда Re 0 при любом шаге h . Размер шага при этом должен выбираться лишь из соображений точности интегрирования.

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

Устойчивость формулы трапеций. Рассмотрим, в заключение,

формулу трапеций

x1 x0 0.5 h ( x0‘ x1‘ ) .

В соответствии с (10.11) подставив значения x0‘ x0 и x1‘ x1 , получим x1 x0 0.5 h ( x0 x1 ) ,

или

x1 x0 (1 0.5 h ) /(1 0.5 h ).

Соответственно для n — го шага, получим

xn [( 2 q ) /( 2 q )] n x0 .

В предельном случае, при n , условие устойчивости имеет вид

( 2 q ) /( 2 q ) 1 ,

или

( 2 u j v ) /( 2 u j v ) 1 .

Преобразуя это условие, получим его в виде 4 u 0 , показывающем, что границей устойчивости в этом случае является мнимая ось, а областью устойчивости — левая полуплоскость рисунок 10.2. Таким образом, формула трапеций устойчива для устойчивых схем, когда Re 0 , при любом шаге интегрирования h .

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

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

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

hi ci x( i )( t0 ) ,

p 1

так называемая ошибка усечения, пропорциональная шагу интегрирования h

Одним из путей обеспечения точности при выбранном шаге h , является использование формулы с наивысшим порядком p . Можно так же

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

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

Следует, наконец, отметить, что полученные нами для простых методов интегрирования границы областей устойчивости справедливы лишь для рассматриваемого дифференциального уравнения вида (10.11). Для других дифференциальных уравнений и других методов интегрирования границы областей будут другими. Исследование границ устойчивости, для простых методов интегрирования и частного вида дифференциального уравнения, предпринято с целью, обозначить основные тенденции и дать сравнительную характеристику методов.

Дадим неформальное объяснение, некоторым часто используемым терминам и понятиям.

Формулы интегрирования, основанные на значениях функции и ее производных в предыдущие моменты времени, подобные прямой формуле Эйлера называют иногда явными. Явные формулы, в силу их слабой устойчивости используются главным образом для предсказания начальных значений при использовании других, неявных формул интегрирования.

Неявными формулами интегрирования, аналогичными обратной формуле Эйлера и формуле трапеций, называются формулы, связывающие значение

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

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

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

Линейные многошаговые формулы, интегрирования, содержащие


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

дифференцирования назад (ФДН). Формулы дифференцирования назад обладают повышенной устойчивостью решений и зачастую используются в качестве основных алгоритмов численного интегрирования дифференциальных уравнений.

Формулу интегрирования называют A — устойчивой, если она дает ограниченное решение тестового дифференциального уравнения x‘ x , для произвольных размеров шага и любого числа шагов, при Re 0 . Обратная формула Эйлера и формула трапеций обладают этим свойством.

Область абсолютной устойчивости, какой-либо формулы интегрирования, это часть плоскости q h , в которой интегрирование

дифференциального уравнения x‘ x , при Re 0 , дает ограниченный результат при любом числе шагов. Так прямая формула Эйлера абсолютно устойчива внутри единичной окружности с единичным радиусом и центром в точке ( 1,0 ). Обратная формула Эйлера и формула трапеций абсолютно

устойчивы во всей левой полуплоскости.

Большинство цепей встречающихся на практике, являются устойчивыми. Для линейной цепи это означает, что действительные части

корней характеристических уравнений отрицательны. В связи с этим, говоря

об интегрировании уравнения x‘ x , имеем в виду решение, при Re 0 . Однако встречаются ситуации, например автоколебательные системы, когда схема в рабочей точке неустойчива и в установившемся режиме в ней наблюдаются периодические колебания. При численном интегрировании дифференциальных уравнений таких цепей следует применять методы, учитывающие эти свойства. Обратная формула Эйлера, например, в этой ситуации не подходит, т.к. дает затухающие решения. В этом случае

предпочтительны, формула трапеций и другие методы интегрировании. Жесткая система дифференциальных уравнений, это система имеющая

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

При интегрировании такой системы по прямой формуле Эйлера для обеспечения устойчивости потребуется очень малый шаг, чтобы q h

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

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

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

10.5 Расчет переходных процессов цепей

Как уже отмечалось, рассмотренные нами формулы численного интегрирования применимы как к линейным, так и нелинейным дифференциальным уравнениям.

Уравнения переменных состояния. Если цепь описана линейной системой дифференциальных уравнений в нормальной форме Коши, относительно переменных состояния

x‘ ( t ) A x( t ) B w( t ),	(10.12)
то вычисление x( t ) не вызывает проблем. Так	вычисление переменных

состояния может быть осуществлено на основании полученных нами выражений в матричной форме (10.7 — 10.9).

Непосредственно метод переменных состояния нами не рассматривался ввиду того, что этот метод в настоящее время используется редко из-за сложности алгоритма формирования системы дифференциальных уравнений в нормальной форме Коши для произвольного вида цепей. Заметим только, что метод переменных состояния достаточно широко освещен в литературе. В качестве переменных состояния выступают независимые токи в индуктивностях и напряжения на емкостях. Для формирования системы дифференциальных уравнений в нормальной форме Коши в методе переменных состояния используются топологические уравнения, описывающие дерево графа через матрицу главных сечений и компонентные уравнения ветвей цепи.

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

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

Система узловых уравнений, при заданных номиналах, имеет вид

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

Источник

Локальная и глобальная ошибки дискретизации при численном интегрировании дифференциальных уравнений

Аппроксимируем производную в момент времени t_k соотношением

При такой аппроксимации уравнение (1) примет вид:

Формула (2) известна как прямой метод Эйлера.

Рис.1 Иллюстрация алгоритмов (а) прямого метода Эйлера, (b) обратного метода Эйлера

Обратный метод Эйлера подобен прямому, но есть одно отличие в аппроксимации для производной

Такая аппроксимация дает формулу обратного метода Эйлера:

Трапецеидальный алгоритм имеет вид:

Локальная ошибка определяется как ошибка на шаге алгоритма:

Для m-шаговых алгоритмов предполагается, что предыдущие m точек x_k-i при i=0. m-1 являются точно заданными, т.е.

Если пренебречь зависимостью a_k от x_k, то на единичном интервале времени ошибка усечения будет:

Порядок метода интегрирования.

Дата добавления: 2015-06-12 ; просмотров: 3180 ; Нарушение авторских прав

· ошибка округления, накопившаяся при выполнении арифметических операций.

При этом доминирующей является, как правило, ошибка дискретизации.

где

;

3.2.6. Процедуры численного моделирования
с автоматическим выбором шага

Первый способ реализует правило Рунге [29, 23], при котором ошибка дискретизации определяется по формуле

где – некоторая константа;

– порядок метода интегрирования.

Для метода Эйлера формула Рунге дает

для метода Рунге – Кутты 4-го порядка

При автоматическом выборе шага наиболее популярным является алгоритм «трех зон» [29]:

где – порядок метода интегрирования.

3.2.7. Особенности численного интегрирования
технических систем

Пусть исследуемая система является линейной и может быть описана матричными уравнениями состояния

Оценкой жесткости системы ОДУ обычно считается число обусловленности матрицы Якоби

где – норма матрицы .

3.3. Моделирование гибридных (событийно-управляемых)
технических систем

Гибридное поведение может быть связано со следующими факторами:

· гибридное поведение, обусловленное изменением состава системы.

Все перечисленные факторы являются типичными для технических систем.

Рис. 3.4. Диаграмма состояний и переходов цифровых часов

Рис. 3.5. Карта состояний и карта поведений обрывающегося маятника

Глава 4
Автоматизированное моделирование
технических объектов

Можно отметить следующие факторы, способствующие внедрению систем автоматизированного моделирования:

4.1. Особенности современных систем
автоматизированного моделирования

Современные профессиональные САМ имеют следующие отличительные черты [31]:

· реализация механизмов продвижения модельного времени, основанных как на принципе , так и на принципе ;

· наличие средств, обеспечивающие формирование виртуальных аналогов измерительно-управляющей аппаратуры.

· визуализация результатов во время эксперимента;

· возможность интерактивного вмешательства в ход моделирования;

· возможность использования 2D и 3D анимации.

4.2. Иерархическое проектирование
и многоуровневое моделирование технических систем

4.4. Архитектура программ
автоматизированного моделирования

4.4.1. Графический интерфейс программ математического
моделирования динамических систем

Рис. 4.1. Модульная структура
системы автоматизированного моделирования

· операторно-структурные схемы, принятые в ТАУ;

· функциональные и принципиальные схемы различных физических устройств;

· кинематические схемы механизмов;

· блок-схемы алгоритмов и другие графические модели.

· преобразование информации о схеме в команды для моделирующей программы (моделятора);

· контроль за процессом моделирования, визуализация результатов моделирования.

П.5. Локальные и глобальные погрешности одношаговых методов решения ДУ

(метода Эйлера и методов Рунге-Кутта 2го, 4го порядка).

Теорема 6.1:

Если локальная погрешность метода , то глобальная .

Методы	Локальная	Глобальная
Эйлер	const*h 2	const*h
Р.–К. 2го порядка по времени	const*h 3	const*h 2
Р.–К. 2го порядка по производной	const*h 3	const*h 2
Р.–К. 4го порядка	const*h 5	const*h 4

Как и при численном интегрировании, порядок метода – степень h в глобальной погрешности.

П.6. Многошаговые методы решения ДУ и СДУ.

В k-шаговом методе используются значения с k предыдущих шагов.

Формулы метода Милна:

(6.9)

П.7. Оценка погрешности решения ДУ и СДУ методом двойного пересчета. Коррекция решения.

Как и при ЧИ, при решении ДУ и СДУ после 2го пересчёта в качестве точного решения выгодно брать не , а .

— для второго порядка

П.7. Краевые задачи для дифференциальных уравнений.

Рассмотрим простейшую краевую задачу для ДУ 2го порядка:

(6.10)

А мы умеем решать:

(6.11),

В (6.11) нам известно , поэтому для решения задачи (6.10) мы будем подбирать в (6.11), с тем, чтобы у(b) = у₁

Метод стрельб

После пристрелки и определения интервала [a,b],

где идёт смена знака, запускаем МПД или МХ.

На практике это выглядит так, как будто мы

решаем уравнение , где возвращает

решение задачи Коши (6.11) в точке b при

П.9. Что делать, если ДУ не может быть разрешено относительно старшей производной?

Решать уравнение относительно старшей производной любым методом(Хорд, МПД, Ньютона).

Замечание:

Типы и классификация ошибок численного интегрирования

1. Неустранимая погрешность исходных данных. На нее повлиять нельзя, но учесть необходимо.
2. Локальная ошибка метода. Ее величина определяется порядком р используемого метода и зависит от шага: е,

„ _ , „ , ного изменения ошибки. Возможно, однако, си-

Рнс.7.1. Суммарная ошибка

Теперь разобьем шаг на два, длины —. и вычислим то же значение,

выполнив два шага интегрирования. Будем считать, что величина R при этом нс изменилась. Получим:

Вычитая уравнения, можно получить

И можно считать, что при интегрировании с постоянным шагом И верно соотношение

Таким образом, схема выбора величины /г на шаге для достижения заданной точности Очакова:

1. при известном х, с шагом /г вычислить т_г+/;

источники:

http://life-prog.ru/2_14935_poryadok-metoda-integrirovaniya.html

http://megalektsii.ru/s20945t2.html

http://bstudy.net/718789/ekonomika/tipy_klassifikatsiya_oshibok_chislennogo_integrirovaniya

Источник

Ошибки усечения в численное интегрирование бывают двух видов:

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

Определения

Предположим, у нас есть непрерывное дифференциальное уравнение

$y '= f (t, y), qquad y (t_ {0}) = y_ {0}, qquad t geq t_ {0}$

и мы хотим вычислить приближение $y_n$ истинного решения $y (t_ {n})$ с дискретными временными шагами $t_ {1}, t_ {2}, ldots, t_ {N}$ . Для простоты предположим, что временные шаги равномерно распределены:

$h = t_ {n} -t _ {{n-1}}, qquad n = 1,2, ldots, N.$

Предположим, мы вычисляем последовательность $y_n$ с одношаговым методом формы

$y_ {n} = y _ {{n-1}} + hA (t _ {{n-1}}, y _ {{n-1}}, h, f).$

Функция $А$ называется функция приращения, и может интерпретироваться как оценка наклона ${ displaystyle { frac {y (t_ {n}) - y (t_ {n-1})} {h}}}$ .

Ошибка локального усечения

В локальная ошибка усечения $тау _ {п}$ ошибка, из-за которой наша функция приращения, $А$ , вызывает во время одной итерации при условии полного знания истинного решения на предыдущей итерации.

Более формально локальная ошибка усечения, $тау _ {п}$ , на шаге $п$ вычисляется из разницы между левой и правой частями уравнения для приращения $y_ {n} приблизительно y _ {{n-1}} + hA (t _ {{n-1}}, y _ {{n-1}}, h, f)$ :

$tau _ {n} = y (t_ {n}) - y (t _ {{n-1}}) - hA (t _ {{n-1}}, y (t _ {{n-1}}), з, е).$ ^[1]^[2]

Численный метод последовательный если локальная ошибка усечения $ой)$ (это означает, что для каждого $varepsilon> 0$ существует $ЧАС$ такой, что $| tau _ {n} | < varepsilon h$ для всех $h <H$ ; увидеть маленькая нотация ). Если функция приращения $А$ непрерывна, то метод непротиворечив тогда и только тогда, когда $A (t, y, 0, f) = f (t, y)$ .^[3]

Кроме того, мы говорим, что численный метод имеет порядок $п$ если для любого достаточно гладкого решения задачи начального значения локальная ошибка усечения равна $О (ч ^ {{р + 1}})$ (это означает, что существуют константы $C$ и $ЧАС$ такой, что $| tau _ {n} | <Ch ^ {{p + 1}}$ для всех $h <H$ ).^[4]

Глобальная ошибка усечения

В глобальная ошибка усечения это накопление локальная ошибка усечения на всех итерациях, предполагая полное знание истинного решения на начальном временном шаге.^{[нужна цитата ]}

Более формально, глобальная ошибка усечения, $e_ {n}$ , вовремя $t_n$ определяется:

${ begin {align} e_ {n} & = y (t_ {n}) - y_ {n} & = y (t_ {n}) - { Big (} y_ {0} + hA (t_ { 0}, y_ {0}, h, f) + hA (t_ {1}, y_ {1}, h, f) + cdots + hA (t _ {{n-1}}, y _ {{n-1 }}, h, f) { Big)}. end {align}}$ ^[5]

Численный метод сходящийся если глобальная ошибка усечения становится равной нулю, когда размер шага становится равным нулю; другими словами, численное решение сходится к точному решению: $lim _ {{h to 0}} max _ {n} | e_ {n} | = 0$ .^[6]

Связь между локальными и глобальными ошибками усечения

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

Глобальная ошибка усечения удовлетворяет рекуррентному соотношению:

$e _ {{n + 1}} = e_ {n} + h { Big (} A (t_ {n}, y (t_ {n}), h, f) -A (t_ {n}, y_ {n) }, h, f) { Big)} + tau _ {{n + 1}}.$

Это сразу следует из определений. Теперь предположим, что функция приращения равна Липшицева непрерывная во втором аргументе, то есть существует постоянная $L$ такой, что для всех $т$ и $г_ {1}$ и $г_ {2}$ , у нас есть:

$| A (t, y_ {1}, h, f) -A (t, y_ {2}, h, f) | leq L | y_ {1} -y_ {2} |.$

Тогда глобальная ошибка удовлетворяет оценке

$| e_ {n} | leq { frac { max _ {j} tau _ {j}} {hL}} left ({ mathrm {e}} ^ {{L (t_ {n} -t_ {0})}} - 1 вправо).$ ^[7]

Из приведенной выше оценки глобальной ошибки следует, что если функция $ж$ в дифференциальном уравнении непрерывна по первому аргументу и липшицева по второму аргументу (условие из Теорема Пикара – Линделёфа ), а функция приращения $А$ непрерывна по всем аргументам и липшицева по второму аргументу, то глобальная ошибка стремится к нулю как размер шага $час$ стремится к нулю (иными словами, численный метод сходится к точному решению).^[8]

Расширение линейных многоступенчатых методов

Теперь рассмотрим линейный многоступенчатый метод, задаваемый формулой

$begin {align} & y_ {n + s} + a_ {s-1} y_ {n + s-1} + a_ {s-2} y_ {n + s-2} + cdots + a_0 y_n & qquad {} = h bigl (b_s f (t_ {n + s}, y_ {n + s}) + b_ {s-1} f (t_ {n + s-1}, y_ {n + s}) -1}) + cdots + b_0 f (t_n, y_n) bigr), end {align}$

Таким образом, следующее значение численного решения вычисляется согласно

${ displaystyle y_ {n + s} = - sum _ {k = 0} ^ {s-1} a_ {k} y_ {n + k} + h sum _ {k = 0} ^ {s} b_ {k} f (t_ {n + k}, y_ {n + k}).}$

Следующая итерация линейного многоступенчатого метода зависит от предыдущего s повторяется. Таким образом, в определении локальной ошибки усечения теперь предполагается, что предыдущие s все итерации соответствуют точному решению:

${ displaystyle tau _ {n} = y (t_ {n + s}) + sum _ {k = 0} ^ {s-1} a_ {k} y (t_ {n + k}) - час сумма _ {k = 0} ^ {s} b_ {k} f (t_ {n + k}, y (t_ {n + k})).}$ ^[9]

Опять же, метод непротиворечив, если $тау _ {п} = о (ч)$ и в нем порядок п если $tau _ {n} = O (h ^ {{p + 1}})$ . Не изменилось и определение глобальной ошибки усечения.

Связь между локальными и глобальными ошибками усечения немного отличается от более простой настройки одношаговых методов. Для линейных многоступенчатых методов существует дополнительная концепция, называемая нулевая стабильность необходим для объяснения связи между локальными и глобальными ошибками усечения. Линейные многоступенчатые методы, удовлетворяющие условию нулевой устойчивости, имеют такое же соотношение между локальными и глобальными ошибками, что и одношаговые методы. Другими словами, если линейный многоступенчатый метод устойчив к нулю и непротиворечив, то он сходится. И если линейный многоступенчатый метод устойчив к нулю и имеет локальную ошибку $tau _ {n} = O (h ^ {{p + 1}})$ , то его глобальная ошибка удовлетворяет $e_ {n} = O (h ^ {p})$ .^[10]

Смотрите также

Порядок точности
Численное интегрирование
Численные обыкновенные дифференциальные уравнения
Ошибка усечения

Заметки

^ Gupta, G.K .; Sacks-Davis, R .; Тишер П. Э. (март 1985 г.). «Обзор последних достижений в решении ODE». Вычислительные опросы. 17 (1): 5–47. CiteSeerX 10.1.1.85.783. Дои:10.1145/4078.4079.
^ Сюли и Майерс 2003, п. 317, звонки $tau _ {n} / ч$ ошибка усечения.
^ Сюли и Майерс 2003, стр. 321 и 322
^ Изерлес 1996, п. 8; Сюли и Майерс 2003, п. 323
^ Сюли и Майерс 2003, п. 317
^ Изерлес 1996, п. 5
^ Сюли и Майерс 2003, п. 318
^ Сюли и Майерс 2003, п. 322
^ Сюли и Майерс 2003, п. 337, использует другое определение, разделив его по существу на час
^ Сюли и Майерс 2003, п. 340

использованная литература

Изерлес, Арье (1996), Первый курс численного анализа дифференциальных уравнений, Издательство Кембриджского университета, ISBN 978-0-521-55655-2.
Сули, Эндре; Майерс, Дэвид (2003), Введение в численный анализ, Издательство Кембриджского университета, ISBN 0521007941.

внешние ссылки

Примечания к ошибкам усечения и методам Рунге-Кутта
Ошибка усечения метода Эйлера

Источник

Локальная и глобальная ошибки дискретизации при численном интегрировании дифференциальных уравнений

Порядок метода интегрирования.

П.5. Локальные и глобальные погрешности одношаговых методов решения ДУ

Типы и классификация ошибок численного интегрирования

Интегрирование функции, заданной таблично

Локальная и глобальная ошибки дискретизации при численном интегрировании дифференциальных уравнений

Порядок метода интегрирования.

П.5. Локальные и глобальные погрешности одношаговых методов решения ДУ

Типы и классификация ошибок численного интегрирования

Определения

Ошибка локального усечения

Глобальная ошибка усечения

Связь между локальными и глобальными ошибками усечения

Расширение линейных многоступенчатых методов

Смотрите также

Заметки

использованная литература

внешние ссылки

Возможно, вам также будет интересно: