Математическое моделирование в технических системах

ЛЕКЦИИ «МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ В ТЕХНИЧЕСКИХ СИСТЕМАХ» 1 ОГЛАВЛЕНИЕ Лекция 1. Математическое моделирование. Форма и принципы представления математических моделей............................................................................................... 3 Лекция 2. Особенности построения математических моделей ................................. 9 Лекция 3. Компьютерное моделирование и вычислительный эксперимент. Решение математических моделей............................................................................. 14 Лекция 4. Компьютерное моделирование и решение линейных многомерных систем ............................................................................................................................ 18 Лекция 5. Моделирование и решение многомерных систем нелинейных уравнений...................................................................................................................... 28 Лекция 6. Компьютерное моделирование при обработке ....................................... 35 опытных данных .......................................................................................................... 35 Лекция 7. Компьютерное моделирование и решение дифференциальных уравнений...................................................................................................................... 53 2 Лекция 1. Математическое моделирование. Форма и принципы представления математических моделей Общие вопросы математического моделирования. Приведена классификация математических моделей ЭВМ прочно вошла в нашу жизнь, и практически нет такой области человеческой деятельности, где не применялась бы ЭВМ. ЭВМ сейчас широко используется в процессе создания и исследования новых машин, новых технологических процессов и поиске их оптимальных вариантов; при решении экономических задач, при решении задач планирования и управления производством на различных уровнях. Создание же крупных объектов в ракетостроении, авиастроении, судостроении, а также проектирование плотин, мостов, и др. вообще невозможно без применения ЭВМ. Для использования ЭВМ при решении прикладных задач, прежде всего прикладная задача должна быть "переведена" на формальный математический язык, т.е. для реального объекта, процесса или системы должна быть построена его математическая модель. Слово "Модель" происходит от латинского modus (копия, образ, очертание). Моделирование - это замещение некоторого объекта А другим объектом Б. Замещаемый объект А называется оригиналом или объектом моделирования, а замещающий Б - моделью. Другими словами, модель - это объект-заменитель объекта-оригинала, обеспечивающий изучение некоторых свойств оригинала. Целью моделирования являются получение, обработка, представление и использование информации об объектах, которые взаимодействуют между собой и внешней средой; а модель здесь выступает как средство познания свойств и закономерности поведения объекта. Моделирование широко используется в различных сферах человеческой деятельности, особенно в сферах проектирования и управления, где особенными являются процессы принятия эффективных решений на основе получаемой информации. Модель всегда строится с определенной целью, которая оказывает влияние на то, какие свойства объективного явления оказываются существенными, а какие нет. Модель представляет собой как бы проекцию объективной реальности под определенным углом зрения. Иногда, в зависимости от целей, можно получить ряд проекций объективной реальности, вступающих в противоречие. Это характерно, 3 как правило, для сложных систем, у которых каждая проекция выделяет существенное для определенной цели из множества несущественного. Теорией моделирования является раздел науки, изучающий способы исследования свойств объектов-оригиналов, на основе замещения их другими объектами-моделями. В основе теории моделирования лежит теория подобия. При моделировании абсолютное подобие не имеет места и лишь стремится к тому, чтобы модель достаточно хорошо отображала исследуемую сторону функционирования объекта. Абсолютное подобие может иметь место лишь при замене одного объекта другим точно таким же. Все модели можно разделить на два класса: 1. вещественные, 2. идеальные. 3. В свою очередь вещественные модели можно разделить на: 4. натурные, 5. физические, 6. математические. Идеальные модели можно разделить на: 1. наглядные, 2. знаковые, 3. математические. Вещественные натурные модели - это реальные объекты, процессы и системы, над которыми выполняются эксперименты научные, технические и производственные. Вещественные физические модели - это макеты, муляжи, воспроизводящие физические свойства оригиналов (кинематические, динамические, гидравлические, тепловые, электрические, световые модели). Вещественные математические это аналоговые, структурные, геометрические, графические, цифровые и кибернетические модели. Идеальные наглядные модели - это схемы, карты, чертежи, графики, графы, аналоги, структурные и геометрические модели. Идеальные знаковые модели - это символы, алфавит, языки программирования, упорядоченная запись, топологическая запись, сетевое представление. Идеальные математические модели - это аналитические, функциональные, имитационные, комбинированные модели. В приведенной классификации некоторые модели имеют двойное толкование (например - аналоговые). Все модели, кроме натурных, можно 4 объединить в один класс мысленных моделей, т.к. они являются продуктом абстрактного мышления человека. Остановимся на одном из наиболее универсальных видов моделирования математическом, ставящим в соответствие моделируемому физическому процессу систему математических соотношений, решение которой позволяет получить ответ на вопрос о поведении объекта без создания физической модели, часто оказывающейся дорогостоящей и неэффективной. Математическое моделирование - это средство изучения реального объекта, процесса или системы путем их замены математической моделью, более удобной для экспериментального исследования с помощью ЭВМ. Математическая модель является приближенным представлением реальных объектов, процессов или систем, выраженным в математических терминах и сохраняющим существенные черты оригинала. Математические модели в количественной форме, с помощью логико-математических конструкций, описывают основные свойства объекта, процесса или системы, его параметры, внутренние и внешние связи. В общем случае математическая модель реального объекта, процесса или системы представляется в виде системы функционалов Фi (X,Y,Z,t)=0, где X - вектор входных переменных, X=[x1,x2,x3, ... , xN]t, Y - вектор выходных переменных, Y=[y1,y2,y3, ... , yN]t, Z - вектор внешних воздействий, Z=[z1,z2,z3, ... , zN]t, t - координата времени. Построение математической модели заключается в определении связей между теми или иными процессами и явлениями, создании математического аппарата, позволяющего выразить количественно и качественно связь между теми или иными процессами и явлениями, между интересующими специалиста физическими величинами, и факторами, влияющими на конечный результат. Обычно их оказывается настолько много, что ввести в модель всю их совокупность не удается. При построении математической модели перед исследованием возникает задача выявить и исключить из рассмотрения факторы, несущественно влияющие на конечный результат (математическая модель обычно включает значительно меньшее число факторов, чем в реальной действительности). На основе данных эксперимента выдвигаются гипотезы о связи между величинами, выражающими конечный результат, и факторами, введенными в математическую модель. Такая связь зачастую выражается системами дифференциальных уравнений в частных производных (например, в задачах 5 механики твердого тела, жидкости и газа, теории фильтрации, теплопроводности, теории электростатического и электродинамического полей). Конечной целью этого этапа является формулирование математической задачи, решение которой с необходимой точностью выражает результаты, интересующие специалиста. Форма и принципы представления математической модели зависит от многих факторов. По принципам построения математические модели разделяют на: 1. аналитические; 2. имитационные. В аналитических моделях процессы функционирования реальных объектов, процессов или систем записываются в виде явных функциональных зависимостей. Аналитическая модель разделяется на типы в зависимости от математической проблемы: 1. уравнения (алгебраические, трансцендентные, дифференциальные, интегральные), 2. аппроксимационные задачи (интерполяция, экстраполяция, численное интегрирование и дифференцирование), 3. задачи оптимизации, 4. стохастические проблемы. Однако по мере усложнения объекта моделирования построение аналитической модели превращается в трудноразрешимую проблему. Тогда исследователь вынужден использовать имитационное моделирование. В имитационном моделировании функционирование объектов, процессов или систем описывается набором алгоритмов. Алгоритмы имитируют реальные элементарные явления, составляющие процесс или систему с сохранением их логической структуры и последовательности протекания во времени. Имитационное моделирование позволяет по исходным данным получить сведения о состояниях процесса или системы в определенные моменты времени, однако прогнозирование поведения объектов, процессов или систем здесь затруднительно. Можно сказать, что имитационные модели - это проводимые на ЭВМ вычислительные эксперименты с математическими моделями, имитирующими поведение реальных объектов, процессов или систем. В зависимости от характера исследуемых реальных процессов и систем математические модели могут быть: 1. детерминированные, 2. стохастические. 6 В детерминированных моделях предполагается отсутствие всяких случайных воздействий, элементы модели (переменные, математические связи) достаточно точно установленные, поведение системы можно точно определить. При построении детерминированных моделей чаще всего используются алгебраические уравнения, интегральные уравнения, матричная алгебра. Стохастическая модель учитывает случайный характер процессов в исследуемых объектах и системах, который описывается методами теории вероятности и математической статистики. По виду входной информации модели разделяются на: 1. непрерывные, 2. дискретные. Если информация и параметры являются непрерывными, а математические связи устойчивы, то модель - непрерывная. И наоборот, если информация и параметры - дискретны, а связи неустойчивы, то и математическая модель дискретная. По поведению моделей во времени они разделяются на: 1. статические, 2. динамические. Статические модели описывают поведение объекта, процесса или системы в какой-либо момент времени. Динамические модели отражают поведение объекта, процесса или системы во времени. По степени соответствия между математической моделью и реальным объектом, процессом или системой математические модели разделяют на: 1. изоморфные (одинаковые по форме), 2. гомоморфные (разные по форме). Модель называется изоморфной, если между нею и реальным объектом, процессом или системой существует полное поэлементное соответствие. Гомоморфной - если существует соответствие лишь между наиболее значительными составными частями объекта и модели. В дальнейшем для краткого определения вида математической модели в приведенной классификации будем пользоваться следующими обозначениями: Первая буква: Д - детерминированная, С - стохастическая. Вторая буква: Н - непрерывная, Д - дискретная. 7 Третья буква: А - аналитическая, И - имитационная. Согласно этим обозначениям, описанная в лекции 2, модель кривошипношатунного механизма (Рис. 2.1.) обозначается как модель вида ДНА (детерминированная, непрерывная, аналитическая), так как: 1. Отсутствует (точнее не учитывается) влияние случайных процессов, т.е. модель детерминированная (Д). 2. Информация и параметры - непрерывные, т.е. модель - непрерывная (Н), 3. Функционирование модели кривошипно-шатунного механизма описано в виде нелинейных трансцендентных уравнений, т.е. модель - аналитическая (А) 8 Лекция 2. Особенности построения математических моделей Процесс построения математической модели. Словесный алгоритм процесса. Для использования ЭВМ при решении прикладных задач прежде всего прикладная задача должна быть "переведена" на формальный математический язык, т.е. для реального объекта, процесса или системы должна быть построена его математическая модель. Математические модели в количественной форме, с помощью логикоматематических конструкций, описывают основные свойства объекта, процесса или системы, его параметры, внутренние и внешние связи. Для построения математической модели необходимо: 1. тщательно проанализировать реальный объект или процесс; 2. выделить его наиболее существенные черты и свойства; 3. определить переменные, т.е. параметры, значения которых влияют на основные черты и свойства объекта; 4. описать зависимость основных свойств объекта, процесса или системы от значения переменных с помощью логико-математических соотношений (уравнения, равенства, неравенства, логико-математические конструкций); 5. выделить внутренние связи объекта, процесса или системы с помощью ограничений, уравнений, равенств, неравенств, логико-математических конструкций; 6. определить внешние связи и описать их с помощью ограничений, уравнений, равенств, неравенств, логико-математических конструкций. Математическое моделирование, кроме исследования объекта, процесса или системы и составления их математического описания, также включает: 1. построение алгоритма, моделирующего поведение объекта, процесса или системы; 2. проверка адекватности модели и объекта, процесса или системы на основе вычислительного и натурного эксперимента; 3. корректировка модели; 4. использование модели. Математическое описание исследуемых процессов и систем зависит от: 9 1. природы реального процесса или системы и составляется на основе законов физики, химии, механики, термодинамики, гидродинамики, электротехники, теории пластичности, теории упругости и т.д. 2. требуемой достоверности и точности изучения и исследования реальных процессов и систем. На этапе выбора математической модели устанавливаются: линейность и нелинейность объекта, процесса или системы, динамичность или статичность, стационарность или нестационарность, а также степень детерминированности исследуемого объекта или процесса. При математическом моделировании сознательно отвлекаются от конкретной физической природы объектов, процессов или систем и, в основном, сосредотачиваются на изучении количественных зависимостей между величинами, описывающими эти процессы. Математическая модель никогда не бывает полностью тождественна рассматриваемому объекту, процессу или системе. Основанная на упрощении, идеализации, она является приближенным описанием объекта. Поэтому результаты, полученные при анализе модели, носят приближенный характер. Их точность определяется степенью адекватности (соответствия) модели и объекта. Построение математической модели обычно начинается с построения и анализа простейшей, наиболее грубой математической модели рассматриваемого объекта, процесса или системы. В дальнейшем, в случае необходимости, модель уточняется, делается ее соответствие объекту более полным. Возьмем простой пример. Нужно определить площадь поверхности письменного стола. Обычно для этого измеряют его длину и ширину, а затем перемножают полученные числа. Такая элементарная процедура фактически обозначает следующее: реальный объект (поверхность стола) заменяется абстрактной математической моделью – прямоугольником. Прямоугольнику приписываются размеры, полученные в результате измерения длины и ширины поверхности стола, и площадь такого прямоугольника приближенно принимается за искомую площадь стола. Однако модель прямоугольника для письменного стола – это простейшая, наиболее грубая модель. При более серьезном подходе к задаче прежде, чем воспользоваться для определения площади стола моделью прямоугольника, эту модель нужно проверить. Проверки можно осуществить следующим образом: измерить длины противоположных сторон стола, а также длины его диагоналей и сравнить их между собой. Если, с требуемой степенью точности, длины противоположных сторон и длины диагоналей попарно равны между собой, то поверхность стола действительно можно рассматривать как прямоугольник. В 10 противном случае модель прямоугольника придется отвергнуть и заменить моделью четырехугольника общего вида. При более высоком требовании к точности может возникнуть необходимость пойти в уточнении модели еще дальше, например, учесть закругления углов стола. С помощью этого простого примера было показано, что математическая модель не определяется однозначно исследуемым объектом, процессом или системой. Для одного и того же стола мы можем принять либо модель прямоугольника, либо более сложную модель четырехугольника общего вида, либо четырехугольника с закругленными углами. Выбор той или иной модели определяется требованием точности. С повышением точности модель приходится усложнять, учитывая новые и новые особенности изучаемого объекта, процесса или системы. Рассмотрим другой пример: исследование движения кривошипно-шатунного механизма (Рис. 2.1). Рис. 2.1. Для кинематического анализа этого механизма, прежде всего, необходимо построить его кинематическую модель. Для этого: 1. Заменяем механизм его кинематической схемой, где все звенья заменены жесткими элементами и связями; 2. Пользуясь этой схемой, мы выводим уравнение движения механизма; 3. Дифференцируя последнее, получаем уравнения скоростей и ускорения, которые представляют собой дифференциальные уравнения 1-го и 2-го порядка. Запишем эти уравнения: где С0 – крайнее правое положение ползуна С: r – радиус кривошипа AB; l – длина шатуна BC; – угол поворота кривошипа; 11 Полученные трансцендентные уравнения представляют математическую модель движения плоского аксиального кривошипно-шатунного механизма, основанную на следующих упрощающих предположениях: 1. нас не интересовали конструктивные формы и расположение масс, входящих в механизм тел, и все тела механизма мы заменили отрезками прямых. На самом деле, все звенья механизма имеют массу и довольно сложную форму. Например, шатун – это сложное сборное соединение, форма и размеры которого, конечно, будут влиять на движение механизма; 2. при построении математической модели движения рассматриваемого механизма мы также не учитывали упругость входящих в механизм тел, т.е. все звенья рассматривали как абстрактные абсолютно жесткие тела. В действительности же, все входящие в механизм тела – упругие тела. Они при движении механизма будут как-то деформироваться, в них могут даже возникнуть упругие колебания. Это все, конечно, также будет влиять на движение механизма; 3. мы не учитывали погрешность изготовления звеньев, зазоры в кинематических парах A, B, C и т.д. Таким образом, важно еще раз подчеркнуть, что, чем выше требования к точности результатов решения задачи, тем больше необходимость учитывать при построении математической модели особенности изучаемого объекта, процесса или системы. Однако, здесь важно во время остановиться, так как сложная математическая модель может превратиться в трудно разрешимую задачу. Наиболее просто строится модель, когда хорошо известны законы, определяющие поведение и свойства объекта, процесса или системы, и имеется большой практический опыт их применения. Более сложная ситуация возникает тогда, когда наши знания об изучаемом объекте, процессе или системе недостаточны. В этом случае при построении математической модели приходится делать дополнительные предположения, которые носят характер гипотез, такая модель называется гипотетической. Выводы, полученные в результате исследования такой гипотетической модели, носят условный характер. Для проверки выводов необходимо сопоставить результаты исследования модели на ЭВМ с результатами натурного эксперимента. Таким образом, вопрос применимости некоторой математической модели к изучению рассматриваемого объекта, процесса или системы не является математическим вопросом и не может быть решен математическими методами. Основным критерием истинности является эксперимент, практика в самом широком смысле этого слова. 12 Построение математической модели в прикладных задачах – один из наиболее сложных и ответственных этапов работы. Опыт показывает, что во многих случаях правильно выбрать модель – значит решить проблему более, чем наполовину. Трудность данного этапа состоит в том, что он требует соединения математических и специальных знаний. Поэтому очень важно, чтобы при решении прикладных задач математики обладали специальными знаниями об объекте, а их партнеры, специалисты, – определенной математической культурой, опытом исследования в своей области, знанием ЭВМ и программирования. 13 Лекция 3. Компьютерное моделирование и вычислительный эксперимент. Решение математических моделей Компьютерное моделирование. Методы решения математических задач. Компьютерное моделирование как новый метод научных исследований основывается на: 1. построении математических моделей для описания изучаемых процессов; 2. использовании новейших вычислительных машин, обладающих высоким быстродействием (миллионы операций в секунду) и способных вести диалог с человеком. Суть компьютерного моделирования состоит в следующем: на основе математической модели с помощью ЭВМ проводится серия вычислительных экспериментов, т.е. исследуются свойства объектов или процессов, находятся их оптимальные параметры и режимы работы, уточняется модель. Например, располагая уравнением, описывающим протекание того или иного процесса, можно изменяя его коэффициенты, начальные и граничные условия, исследовать, как при этом будет вести себя объект. Более того, можно спрогнозировать поведение объекта в различных условиях. Вычислительный эксперимент позволяет заменить дорогостоящий натурный эксперимент расчетами на ЭВМ. Он позволяет в короткие сроки и без значительных материальных затрат осуществить исследование большого числа вариантов проектируемого объекта или процесса для различных режимов его эксплуатации, что значительно сокращает сроки разработки сложных систем и их внедрение в производство. Компьютерное моделирование и вычислительный эксперимент как новый метод научного исследования заставляет совершенствовать математический аппарат, используемый при построении математических моделей, позволяет, используя математические методы, уточнять, усложнять математические модели. Наиболее перспективным для проведения вычислительного эксперимента является его использование для решения крупных научно-технических и социальноэкономических проблем современности (проектирование реакторов для атомных электростанций, проектирование плотин и гидроэлектростанций, магнитогидродинамических преобразователей энергии, и в области экономики – составление сбалансированного плана для отрасли, региона, для страны и др.). 14 В некоторых процессах, где натурный эксперимент опасен для жизни и здоровья людей, вычислительный эксперимент является единственно возможным (термоядерный синтез, освоение космического пространства, проектирование и исследование химических и других производств). Для проверки адекватности математической модели и реального объекта, процесса или системы результаты исследований на ЭВМ сравниваются с результатами эксперимента на опытном натурном образце. Результаты проверки используются для корректировки математической модели или решается вопрос о применимости построенной математической модели к проектированию либо исследованию заданных объектов, процессов или систем. В заключение подчеркнем еще раз, что компьютерное моделирование и вычислительный эксперимент позволяют свести исследование "нематематического" объекта к решению математической задачи. Этим самым открывается возможность использования для его изучения хорошо разработанного математического аппарата в сочетании с мощной вычислительной техникой. На этом основано применение математики и ЭВМ для познания законов реального мира и их использования на практике. В задачах проектирования или исследования поведения реальных объектов, процессов или систем математические модели, как правило, нелинейны, т.к. они должны отражать реальные физические нелинейные процессы, протекающие в них. При этом параметры (переменные) этих процессов связаны между собой физическими нелинейными законами. Поэтому в задачах проектирования или исследования поведения реальных объектов, процессов или систем чаще всего используются математические модели типа ДНА. Согласно классификации приведенной в лекции 1: Д – модель детерминированная, отсутствует (точнее не учитывается) влияние случайных процессов. Н – модель непрерывная, информация и параметры непрерывны. А – модель аналитическая, функционирование модели описывается в виде уравнений (линейных, нелинейных, систем уравнений, дифференциальных и интегральных уравнений). Итак, мы построили математическую модель рассматриваемого объекта, процесса или системы, т.е. представили прикладную задачу как математическую. После этого наступает второй этап решения прикладной задачи – поиск или разработка метода решения сформулированной математической задачи. Метод должен быть удобным для его реализации на ЭВМ, обеспечивать необходимое качество решения. 15 Все методы решения математических задач можно разделить на 2 группы: 1. точные методы решения задач; 2. численные методы решения задач. В точных методах решения математических задач ответ удается получить в виде формул. Например, вычисление корней квадратного уравнения: или, например, вычисление производных функций: или вычисление определенного интеграла: Однако, подставляя числа в формулу в виде конечных десятичных дробей, мы все равно получаем приближенные значения результата. Для большинства задач, встречающихся на практике, точные методы решения или неизвестны, или дают очень громоздкие формулы. Однако, они не всегда являются необходимыми. Прикладную задачу можно считать практически решенной, если мы сумеем ее решить с нужной степенью точности. Для решения таких задач разработаны численные методы, в которых решение сложных математических задач сводится к последовательному выполнению большого числа простых арифметических операций. Непосредственная разработка численных методов относится к вычислительной математике. Примером численного метода является метод прямоугольников для приближенного интегрирования, не требующий вычисления первообразной для подынтегральной функции. Вместо интеграла вычисляется конечная квадратурная сумма: где x1=a – нижний предел интегрирования; xn+1=b – верхний предел интегрирования; n – число отрезков, на которые разбит интервал интегрирования (a,b) ; 16 – длина элементарного отрезка; f(xi) – значение подынтегральной функции на концах элементарных отрезков интегрирования. Чем больше число отрезков n, на которые разбит интервал интегрирования, тем ближе приближенное решение к истинному, т.е. тем точнее результат. Таким образом, в прикладных задачах и при применении точных методов решения, и при применении численных методов решения результаты вычислений носят приближенный характер. Важно только добиться того, чтобы ошибки укладывались в рамки требуемой точности. Численные методы решения математических задач известны давно, еще до появления ЭВМ, но ими пользовались редко и только в сравнительно простых случаях в силу чрезвычайной трудоемкости вычислений. Широкое применение численных методов стало возможным благодаря ЭВМ. 17 Лекция 4. Компьютерное моделирование и решение линейных многомерных систем Методы и алгоритм решения систем линейных уравнений. Метод Гаусса При моделировании инженерно-технических задач, таких как задачи нахождения неизвестных параметров при рассмотрении условий равновесия системы твердых тел, задач определения оптимального размещения оборудования, оптимального плана производства, оптимального плана перевозок грузов (транспортная задача), технических задач распределения кадров и др., , может быть положена гипотеза линейного представления реального мира. Математические модели таких задач представляются линейными уравнениями. Если задача многомерна, то ее математическая модель представляется системой линейных уравнений. Так, например, к решению систем линейных уравнений сводятся задачи расчета прочности методом конечных элементов. Линейные математические модели также используются в нелинейных системах при условии, если эта нелинейная система условно линеаризирована. В общем виде система линейных уравнений имеет вид: где aij - коэффициенты при неизвестных системы, bi - свободные члены, xj - неизвестные системы, - номер строки, - номер столбца, n - порядок системы. В матричной форме система линейных уравнений имеет вид: где 18 Численные методы решения систем линейных уравнений (СЛУ) можно разделить на две группы: 1. точные или прямые методы, 2. приближенные методы. Приближенные методы реализуют на ЭВМ нахождение корней с заданной точностью и являются итерационными методами. Точные методы позволяют получить решение системы за конечное число итераций. К точным методам относятся: • правило Крамера, • метод Гаусса, • метод прогонки. Решение систем линейных уравнений методом Гаусса Метод Гаусса является точным методом. Он позволяет получить решение системы за конечное число арифметических действий. В основе метода лежит идея последовательного исключения неизвестных. Метод состоит из двух этапов. На первом этапе (прямой ход) система при помощи последовательного исключения неизвестных приводится к треугольному виду. На втором этапе (обратный ход) из системы треугольного вида последовательно, в обратном порядке, начиная c n-го уравнения, находятся неизвестные системы. В качестве примера возьмем систему 4 порядка. 19 (4.1) Прямой ход. На первом шаге прямого хода (к=1) находим x1 из первого уравнения системы (4.1). - ведущий элемент первой строки. Если , то (4.2) Обозначим: (4.3) Подставляя (4.3) в (4.2), получим (4.4) где Подставляем (4.4) во 2, 3 и 4 уравнение системы (4.1), получим: Обозначив коэффициенты при неизвестных полученной системы через а свободные члены через , перепишем полученную систему: (4.5) где 20 Таким образом, в результате выполнения первого шага прямого хода исходная система (4.1) n-го порядка преобразована к совокупности уравнения (4.4) и системы линейных уравнений (4.5), порядок которой равен n-1. На втором шаге прямого хода (к=2) из первого уравнения системы (4.5) находим x2. -ведущий элемент первой строки системы (4.5). Если , то из первого уравнения системы (4.5) имеем: (4.6) где Подставив выражение (4.6) во второе и третье уравнения системы (4.5), получим новую систему линейных уравнений, порядок которой равен n-2. (4.7) где Таким образом, в результате выполнения второго шага прямого хода исходная система (4.1) преобразована к совокупности уравнений (4.4), (4.6) и системы линейных уравнений (4.7),порядок которой равен n-2. На третьем шаге прямого хода (к=3) из системы (4.7) находим x3. - ведущий элемент системы (4.7). Если , то из первого уравнения системы (4.7) имеем: (4.8) где Подставив выражение (4.8) для x3 во второе уравнение системы (4.7) получим: 21 (9.9) где На последнем шаге прямого хода, если то из уравнения (4.9) имеем: (4.10) где (4.11) В результате выполнения всех шагов прямого хода исходная система (4.1) приводится к системе треугольного вида, полученной объединением уравнений (4.4), (4.6), (4.8), 4.10): (4.12) При построении алгоритма прямого хода вычисление организуем в цикле по шагам, т.е. . Последний n-й шаг прямого хода выведем из цикла т.к. здесь реализуется только одно вычисление (4.13) В процессе выполнения всех шагов прямого хода все преобразования коэффициентов и свободных членов проводим по полученным ранее рекуррентным формулам: (4.14) где – номер шага прямого хода, - номер уравнения систем (4.5), (4.7) 22 В процессе обратного хода из системы (4.12) неизвестные находятся в обратном порядке. Значение корня х4 находят из последнего уравнения системы (4.12). Далее х4 используется для отыскания корня х3 из 3-го уравнения, далее х3 и х4 используются отыскания х2 из 2-го уравнения системы (4.12), и, наконец, х2, х3 и х4 используются для отыскания х1 из 1-го уравнения системы (4.12). Все вычисления обратного хода проводим в цикле по i, где по рекуррентным формулам: xi= bi. Рассмотренный выше простейший вариант метода Гаусса, называемый схемой единственного деления, обладает следующим недостатком: если ведущий элемент akk какой-либо строки окажется равным нулю, то этот метод формально непригоден, хотя система может иметь единственное решение. Из этих соображений в схеме алгоритма добавлен поиск ненулевого ведущего элемента. На рисунке 4.1 представлена укрупнённая схема алгоритма (блок-схема) метода Гаусса. 23 Рис.4.1. Укрупнённая схема алгоритма (блок-схема) метода Гаусса Блок 2. С помощью двух вложенных циклов с управляющими переменными i=1,n и j=1,k организуем ввод коэффициентов ai,j и свободных членов bi исходной системы. Для того, чтобы в дальнейшем можно было выполнить в блоке 9 проверку результата, в алгоритме предусмотрено сохранение значений ai,j и bi исходной системы с помощью переприсвоений: cij=aij и di=bi Блок 3. Организуем цикл по k, внутри которого производится вычисление по всем шагам прямого хода. Последний п-й шаг прямого хода выводим из цикла. Блок 4. На каждом шаге прямого хода выполняем поиск ненулевого ведущего элемента. 24 Рис. 4.1. Схема алгоритма (блок-схема) метода Гаусса (продолжение) Поиск ненулевого ведущего элемента ведётся в следующем порядке: а) На каждом k-ом шаге прямого хода ведущий элемент каждой строки сравнивается с нулём; б) Если в k-ой строке имеется нулевой ведущий элемент, то в k-ом столбце в цикле осуществляется поиск ненулевого элемента. в) Если в какой-то строке kn такой ненулевой элемент найден, то строки kn и k поэлементно, в цикле по k1=(k+1),n, меняем местами. Для перестановки элементов используется рабочая переменная R. г) Если ненулевой ведущий элемент не найден, то коду ошибки kо присваиваем значение 1 и расчёт прекращается. Блок 5 - шаг прямого хода. На каждом шаге прямого хода проводим исключение неизвестных путём преобразования коэффициентов и свободных членов системы по полученным ранее рекуррентным формулам. Блок 6. В этом блоке выведем из цикла по k последний шаг прямого хода, т.к. на этом шаге не нужны преобразования коэффициентов и свободных членов, а реализуется только одно вычисление 25 xn=bn/an,n Блок 7 - обратный ход. В процессе обратного хода метода Гаусса из системы треугольного вида последовательно в обратном порядке в цикле по i=(n1),1,-1 находим неизвестные системы по рекуррентной формуле bi= bi - xj.ai,j , i=(n-1),1, j=(n+1),n. При этом в цикле по j=(i+1),n использован приём последовательного вычитания xj.ai,j из bi,после чего вводится переприсвоение bi =хi. Рис. 4.1. Схема алгоритма (блок-схема) метода Гаусса (продолжение) Блок 9 - проверка результата. В этом блоке подставляя значения полученных неизвестных в исходную систему и используя сохранённые значения коэффициентов системы ci,j и свободных членов di, проводим проверку решения задачи по формуле Если корни системы найдены, то Fi – это число, близкое к нулю. Блок 9 в алгоритме метода Гаусса рекомендуется использовать только в процессе отладки метода. В дальнейшем, при использовании метода Гаусса при решении различных прикладных задач, особенно в тех случаях, когда метод Гаусса используется внутри другого метода, блок 9 можно опустить, а в блоке 2 при вводе данных исходные значения коэффициентов системы и её свободных членов можно не сохранять. 26 27 Лекция 5. Моделирование и решение многомерных систем нелинейных уравнений В лекции рассматриваются методы решения систем нелинейных уравнений. В задачах проектирования и исследования поведения реальных объектов, процессов и систем (ОПС) математические модели должны отображать реальные физические нелинейные процессы. При этом эти процессы зависят, как правило, от многих переменных. В результате математические системами нелинейных уравнений. модели реальных ОПС описываются Решение систем нелинейных уравнений Дана система нелинейных уравнений (5.1) или Необходимо решить эту систему, т.е. найти вектор удовлетворяющий системе (5.1) с точностью . , Вектор определяет точку в n-мерном Евклидовом пространстве, т.е. этому пространству и удовлетворяет всем уравнениям системы (5.1). В отличие от систем линейных уравнений для систем нелинейных уравнений неизвестны прямые методы решения. При решении систем нелинейных уравнений используются итерационные методы. Эффективность всех итерационных методов зависит от выбора начального приближения (начальной точки), т.е. вектора . Область, в которой начальное приближение сходится к искомому решению, называется областью сходимости G. Если начальное приближение лежит за пределами G, то решение системы получить не удается. Выбор начальной точки специалиста. во многом определяется интуицией и опытом Метод простых итераций Для применения этого метода исходная система (5.1) должна быть преобразована к виду 28 (5.2) или Далее, выбрав начальное приближение и используя систему (10.2), строим итерационный процесс поиска по схеме: т.е. на каждом k-ом шаге поиска вектор переменных находим, используя значения переменных, полученных на шаге (k-1). Итерационный процесс поиска прекращается как только выполнится условие (5.3) При этом условие (5.3) должно выполняться одновременно по всем переменным. Метод простых итераций используется для решения таких систем линейных уравнений, в которых выполняется условие сходимости итерационного процесса поиска, а именно: (5.4) т.е. сумма абсолютных величин частных производных всех преобразованных уравнений системы (5.2) по j-ой переменной меньше единицы. На рисунке 5.1 представлена схема алгоритма решения систем нелинейных уравнений методом простых итераций. Рис. 5.1. Схема алгоритма метода простых итераций 29 Рассмотрим пример. Дана система нелинейных уравнений: Необходимо определить область сходимости системы, выбрать начальную точку и найти одно из решений системы. 1. 2. Строим графики уравнений: Рис. 5.2. Преобразуем систему для решения методом итераций Проверяем условие сходимости (5.4). Для заданной системы оно имеет вид: Находим: В результате условие (5.4) будет иметь вид: Определяем область сходимости G. 30 Граница области сходимости определится при решении системы, Отсюда х1=0,5 ; . В результате область сходимости определится при и . На графике уравнений строим область сходимости G: Рис. 5.3. , принадлежащую области Выбираем начальную точку сходимости G. Используя выбранную начальную точку заданную систему нелинейных уравнений. решаем Решение систем нелинейных уравнений методом Ньютона Дана система нелинейных уравнений (5.5) или Необходимо решить эту систему, т.е. найти вектор , удовлетворяющий системе (5.5) с точностью . Метод Ньютона наиболее распространенный метод решения систем нелинейных уравнений. Он обеспечивает более быструю сходимость по сравнению с методом простых итераций. 31 В основе метода Ньютона лежит идея линеаризации всех нелинейных уравнений системы (55). Сообщим всей системе (5.5) малые приращения hj и разложим каждое уравнение системы (5.5) в ряд Тейлора: (5.6) где hj - приращение по каждой xj; Ri - остаточные нелинейные члены второго и более высоких порядков каждого ряда Тейлора. Если приращения hj таковы, что переменные xj принимают значения близкие к корню, то будем считать, что левые части уравнений системы (5.6) обращаются в нули. Тогда отбросив Ri сведем задачу решения системы нелинейных уравнений (10.5) к решению системы линейных уравнений, в которой неизвестными являются приращения hj, (5.7) Система (5.7) – система линейных уравнений с неизвестными hj, Запишем (5.7) в матричной форме . где (5.8) 32 (10.9) Матрица А, составленая из частных производных ; называется матрицей Якоби или Якобианом. Метод Ньютона состоит из двух этапов: На первом этапе реализации метода Ньютона необходимо построить систему (5.3). На втором этапе, начиная с начальной точки , необходимо решать систему (5.7) на каждом шаге итерационного процесса поиска методом Гаусса. Найденные значения приращений hj используются как поправки к решению, полученному на предыдущем шаге поиска, т.е. (5.10) или Итерационный процесс прекращается, как только выполнится условие (5.11) по всем приращениям одновременно. Определение матрицы Якоби В методе Ньютона на каждом шаге итерационного процесса поиска необходимо формировать матрицу Якоби, при этом каждый элемент матрицы можно определить: 1. аналитически, как частную производную , 2. методом численного дифференцирования, как отношение приращения функции к приращению аргумента, т.е. В результате частная производная определится как а частная производная по первой координате х1 по координате хj определится как 33 где . Метод Ньютона имеет преимущества по сравнению с другими методами. Но для метода Ньютона так же существует проблема сходимости, с увеличением числа неизвестных область сходимости уменьшается, а в случае больших систем, сходимость обеспечивается если начальная точка близка к искомому решению. На рисунке 5.4 представлена укрупнённая схема алгоритма (блок-схема) метода Ньютона. На рисунках 5.5 и 5.6 представлены схемы алгоритмов метода Ньютона с различными способами определения матрицы Якоби. Рис. 5.4. Блок-схема алгоритма метода Ньютона 34 Лекция 6. Компьютерное моделирование при обработке опытных данных В лекции рассматриваются интерполяции опытных данных методы решения задач аппроксимации и Любому специалисту в своей практической деятельности приходится изучать зависимости между различными параметрами исследуемых объектов, процессов и систем. Например: зависимость числа оборотов двигателя от нагрузки, т.е. n=f(Мкр.) ; зависимость силы резания при обработке детали на металлорежущем станке от глубины резания, т.е. P=f(t), и т.д. Из всех способов задания зависимостей наиболее удобным является аналитический способ задания зависимости в виде функции n=f(Мкр.), P=f(t), y=f(t). Однако на практике специалист чаще всего получает зависимости между исследуемыми параметрами экспериментально. В этом случае ставится натурный эксперимент, изменяются значения параметров на входе системы, измеряются значения параметров на выходе системы. Результаты измерений заносятся в таблицу. Таким образом, в результате проведения натурного эксперимента получаем зависимости между исследуемыми параметрами в виде таблицы, т.е. получаем, так называемую, табличную функцию. Далее с этой табличной функцией необходимо вести научноисследовательские расчеты. Например, необходимо проинтегрировать или продифференцировать табличную функцию и т.д. Рассмотрим две задачи по обработке опытных данных: 1. задачу интерполирования, 2. задачу аппроксимации. Интерполирование функций Дана табличная функция, т.е. дана таблица, в которой для некоторых дискретных значений аргумента xi, расположенных в порядке возрастания, заданы соответствующие значения функции уi: i x y x0 y0 1 x1 y1 35 2 x2 y2 ... ... ... i xi yi ... ... ... n xn yn (6.1) Точки с координатами (xi, yi) называются узловыми точками или узлами. Количество узлов в табличной функции равно N=n+1. На графике табличная функция представляется в виде совокупности узловых точек (рис. 6.1). Рис. 6.1. Длина участка [x0, xn] равна (xn - x0). В расчетной практике инженера часто возникают задачи найти значение функции для аргументов, которые отсутствуют в таблице. Такие задачи называются задачами интерполирования или экстраполирования. Задача интерполирования функции (или задача интерполяции) состоит в том, чтобы найти значения yk табличной функции в любой промежуточной точке хк, расположенной внутри интервала [x0, xn], т.е. и 36 Задача экстраполирования функции (или задача экстраполяции) состоит в том, чтобы найти значения yl табличной функции в точке хl, которая не входит в интервал [x0, xn], т.е. Такую задачу часто называют задачей прогноза. Обе эти задачи решаются при помощи нахождения аналитического выражения некоторой вспомогательной функции F(x), которая приближала бы заданную табличную функцию, т.е. в узловых точках принимала бы значение табличных функций Для определенности задачи искомую функцию F(x) будем искать из класса алгебраических многочленов: (6.2) Этот многочлен должен пройти через все узловые точки, т.е. (6.3) Поэтому степень многочлена n зависит от количества узловых точек N и равна количеству узловых точек минус один, т.е. n=N-1. Многочлен вида (6.2), который проходит через все узловые точки табличной функции называется интерполяционным многочленом. Интерполирование с помощью алгебраических многочленов называется параболическим интерполированием. Таким образом, для решения задачи интерполирования прежде всего необходимо решить задачу, которую можно сформулировать следующим образом: для функции , заданной таблично, построить интерполяционный многочлен степени n, который проходит через все узловые точки таблицы: где n -степень многочлена, равная количеству узловых точек N минус один,т.е. n=N-1. В результате, в любой другой промежуточной точке хk, расположенной внутри отрезка [x0,xn] выполняется приближенное равенство Pn(xk) = f(xk) = yk. (рис.6.2) 37 Рис. 6.2. Построение интерполяционного многочлена в явном виде Для построения интерполяционного многочлена вида (6.2) необходимо определить его коэффициенты a0, a1, :, an, т.е. ai i=0,1,2,:,n. Количество неизвестных коэффициентов равно n+1=N, где n-степень многочлена (6.2), N-количество узловых точек табличной функции (6.1). Для нахождения коэффициентов, используем интерполяционного многочлена (6.2). На основании свойство (6.3) этого свойства интерполяционный многочлен должен пройти через каждую узловую точку (xi, yi) таблицы (6.1), т.е., (6.4) Подставляя в (6.4) каждую узловую точку таблицы (5.1) получаем систему линейных уравнений: (6.5) Неизвестными системы (6.5) являются a0, a1, a2, :, an т.е. коэффициенты многочлена (6.2). Коэффициенты при неизвестных системы (6.5) легко могут быть определены на основании данных таблицы (6.1). Интерполяция по Лагранжу Интерполяционный многочлен может быть построен при помощи специальных интерполяционных формул Лагранжа, Ньютона, Стерлинга, Бесселя и др. 38 Интерполяционный многочлен по формуле Лагранжа имеет вид: (6.6) Таким образом, интерполяционный многочлен Лагранжа приближает заданную табличную функцию, т.е. Ln(xi) = yi и мы можем использовать его в качестве вспомогательной функции для решения задач интерполирования, т.е. . Чем больше узлов интерполирования на отрезке [x0,xn], тем точнее интерполяционный многочлен приближает заданную табличную функцию (6.1), т.е. тем точнее равенство: Однако с увеличением числа узлов интерполирования возрастает степень интерполяционного многочлена n и в результате значительно возрастает объем вычислительной работы. Поэтому при большом числе узлов необходимо применять ЭВМ. В этом случае удобно находить значения функции в промежуточных точках, не получая многочлен в явном виде. При решении задачи экстраполирования функции с помощью интерполяционного многочлена вычисление значения функции за пределами отрезка [x0,xn] обычно производят не далее, чем на один шаг h, равный наименьшей величине так как за пределами отрезка [x0,xn] погрешности, как правило, увеличиваются. Программирование формулы Лагранжа Свернем формулу Лагранжа (6.6). В результате получим где но при этом обязательно выполнение условия . При построении алгоритма используют конструкцию из двух включенных циклов: 39 Внешним циклом накапливаем сумму . Внутренним циклом накапливаем произведение . Алгоритм (рис.6.3) не предусматривает получение интерполяционного многочлена в явном виде, а сразу решает задачу интерполирования функции в заданной точке, x=D. Обозначения в алгоритме: n - степень интерполяционного многочлена Лагранжа (6.6), равная количеству узловых точек N минус один, т.е. n=N-1. D - значение аргумента в точке, для интерполирования табличной функции (6.1). L - значение многочлена (6.6). Рис. 6.3. Схема алгоритма интерполяции по Лагранжу 40 которой решается задача Интерполяция по Ньютону Дана табличная функция: i xi yi x0 y0 1 x1 y1 2 x2 y2 ... ... ... n xn yn или (6.1) Точки с координатами (xi, yi) называются узловыми точками или узлами. Количество узлов в табличной функции равно N=n+1. Необходимо найти значение этой функции в промежуточной точке, например, x=D, причем . Для решения задачи строим интерполяционный многочлен. Интерполяционный многочлен по формуле Ньютона имеет вид: (6.7) где n - степень многочлена, - разделенные разности 0го, 1-го, 2-го,:., n-го порядка, соответственно. Разделенные разности Значения f(x0), f(x1), : , f(xn), т.е. значения табличной функции в узлах, называются разделенными разностями нулевого порядка (k=0). Отношение называется разделенной разностью первого порядка (k=1) на участке [x0, x1] и равно разности разделенных разностей нулевого порядка на концах участка [x0, x1], разделенной на длину этого участка. Для произвольного участка [xi, xi+1] разделенная разность первого порядка (k=1) равна 41 Отношение называется разделенной разностью второго порядка (k=2) на участке [x0, x2] и равно разности разделенных разностей первого порядка, разделенной на длину участка [x0, x2]. Для произвольного участка [xi, xi+2] разделенная разность второго порядка (k=2) равна Таким образом, разделенная разность k -го порядка на участке [xi, xi+k] может быть определена через разделенные разности (k-1) -го порядка по рекуррентной формуле: 6.8) где n - степень многочлена. Максимальное значение k равно n. Тогда i =0 и разделенная разность n -го порядка на участке [x0,xn] равна , т.е. равна разности разделенных разностей (n-1) -го порядка, разделенной на длину участка [x0,xn]. Разделенные разности являются вполне определенными числами, поэтому выражение (6.7) действительно является алгебраическим многочленом n -й степени. При этом в многочлене (6.7) . все разделенные разности определены для участков [x0, x0+k], Лемма: алгебраический многочлен (6.7), построенный по формулам Ньютона, действительно является интерполяционным многочленом, т.е. значение многочлена в узловых точках равно значению табличной функции Заметим, что решение задачи интерполяции по Ньютону имеет некоторые преимущества по сравнению с решением задачи интерполяции по Лагранжу. Каждое слагаемое интерполяционного многочлена Лагранжа зависит от всех значений табличной функции yi, i=0,1,:n. Поэтому при изменении количества узловых точек N и степени многочлена n (n=N-1) интерполяционный многочлен Лагранжа требуется строить заново. В многочлене Ньютона при изменении 42 количества узловых точек N и степени многочлена n требуется только добавить или отбросить соответствующее число стандартных слагаемых в формуле Ньютона (6.7). Это удобно на практике и ускоряет процесс вычислений. Программирование формулы Ньютона Для построения многочлена Ньютона по формуле (6.7) организуем циклический вычислительный процесс по . При этом на каждом шаге поиска находим разделенные разности k -го порядка. Будем помещать разделенные разности на каждом шаге в массив Y. Тогда рекуррентная формула (11.8) будет иметь вид: (6.9) В формуле Ньютона (6.7) используются разделенные разности k -го порядка, подсчитанные только для участков [x0, x0+k], т.е. разделенные разности k -го порядка для i=0. Обозначим эти разделенные разности k-го порядка как у0. А разделенные разности, подсчитанные для I > 0, используются для расчетов разделенных разностей более высоких порядков. Используя (6.9), свернем формулу (6.7). В результате получим (6.10) где у0 - значение табличной функции (6.1) для x=x0. - разделенная разность k-го порядка для участка [x0, x0+k]. Для вычисления Р удобно использовать рекуррентную формулу P = P(x - xk-1) внутри цикла по k. Схема алгоритма интерполяции по Ньютону представлена на рис.6.4. 43 Рис. 11.4. Схема алгоритма интерполяции по Ньютону Пример интерполяции по Ньютону Дана табличная функция: i xi yi 2 0,693147 1 3 1,098613 2 4 1,386295 3 5 1,609438 Вычислить разделенные разности 1-го, 2-го, 3-го порядков (n=3) и занести их в диагональную таблицу. Разделенные разности первого порядка: 44 Разделенные разности второго порядка: Разделенная разность третьего порядка: Таблица 6.1. Диагональная таблица разделенных разностей Разделенная разность i 0-го пор. 1-го пор. 2-го пор. 3-го пор. 0,693147 0,405466 1,098613 -0,058892 0,287682 0,00887416 1,386295 -0,0322695 0,223143 1,60943 Интерполяционный многочлен Ньютона для заданной табличной функции имеет вид: Далее полученный интерполяционный многочлен Ньютона можно привести к нормальному виду и использовать его для решения задач интерполирования или прогноза. Сплайн-интерполяция Сплайны стали широко использоваться в вычислительной математике сравнительно недавно. В машиностроительном черчении они применяются уже давно, так как сплайны - это лекала или гибкие линейки, деформация которых позволяет провести кривую через заданные точки (xi, уi). Используя теорию изгиба бруса при малых деформациях, можно показать, что сплайн - это группа кубических многочленов, в местах сопряжения которых первая и вторая производные непрерывны. Такие функции называются кубическими сплайнами. Для их построения необходимо задать коэффициенты, 45 которые единственным образом определяют многочлен в промежутке между данными точками. Например, для некоторых функций (рис.6.5) необходимо задать все кубические функции q1(x), q2(x), :qn(x). В наиболее общем случае эти многочлены имеют вид: где kij - коэффициенты, определяемые описанными ранее условиями, количество которых равно 4n. Для определения коэффициентов kij необходимо построить и решить систему порядка 4n. Рис. 6.5. Первые 2n условий требуют, чтобы сплайны соприкасались в заданных точках: Следующие (2п-2) условий требуют, чтобы в местах соприкосновения сплайнов были равны первые и вторые производные: Система алгебраических уравнений имеет решение, если число уравнений соответствует числу неизвестных. Для этого необходимо ввести еще два уравнения. Обычно используются следующие условия: При построении алгоритма метода первые и вторые производные удобно аппроксимировать разделенными разностями соответствующих порядков. 46 Полученный таким образом сплайн называется естественным кубическим сплайном. Найдя коэффициенты сплайна, используют эту кусочно-гладкую полиноминальную функцию для представления данных при интерполяции. Аппроксимация опытных данных В результате проведения натурного эксперимента получена табличная функция: i X Y xo yo 1 x1 y1 2 x2 y2 3 x3 y3 : : : n xn yn Где N-количество узловых точек в таблице, n=N-1. Задача аппроксимации заключается в отыскании аналитической зависимости y=f(x) полученной табличной функции. В настоящее время существует 2 способа аппроксимации опытных данных: Первый способ. Этот способ требует, чтобы аппроксимирующая кривая F(x), аналитический вид которой необходимо найти, проходила через все узловые точки таблицы. Эту задачу можно решить с помощью построения интерполяционного многочлена степени n: (6.12) Однако этот способ аппроксимации опытных данных имеет недостатки: 1. Точность аппроксимации гарантируется в небольшом интервале [x0, xn] при количестве узловых точек не более 7-8. 2. Значения табличной функции в узловых точках должны быть заданы с большой точностью. Известно, что как бы точно не проводился эксперимент, результаты эксперимента содержат погрешности. Дело в том, что на самом деле исследуемая величина зависит не только от одного аргумента Х, но и от других случайных факторов, которые от опыта к опыту колеблются по своим собственным случайным законам. Этим самым обуславливается случайная колеблемость исследуемой функции. 47 В результате аппроксимировать опытные данные с помощью интерполяционного многочлена, который проходил бы через все узловые точки таблицы, не всегда удается. Более того, стремясь пройти через все узловые точки таблицы и увеличивая порядок многочлена, мы тем самым начинаем воспроизводить не только закономерные изменения снимаемой функции, но и ее случайные помехи. Второй способ. На практике нашел применение другой способ аппроксимации опытных данных - сглаживание опытных данных. Сущность этого метода состоит в том, что табличные данные аппроксимируют кривой F(x), которая не обязательно должна пройти через все узловые точки, а должна как бы сгладить все случайные помехи табличной функции. Сглаживание опытных данных методом наименьших квадратов В этом методе при сглаживании опытных данных аппроксимирующей кривую F(x) стремятся провести так, чтобы ее отклонения от табличных данных (уклонения) по всем узловым точкам были минимальными (рис 11.6), т.е. (11.6) Рис. 6.6. Избавимся от знака уклонения. Тогда условие (6.6) будет иметь вид: (6.7) Суть метода наименьших квадратов заключается в следующем: для табличных данных, полученных в результате эксперимента, отыскать аналитическую зависимость F(x), сумма квадратов уклонений которой от табличных данных по всем узловым точкам была бы минимальной, т.е. (6.8) Для определенности задачи искомую функцию F(x) будем выбирать из класса алгебраических многочленов степени m: 48 (6.9) Назовем многочлен аппроксимирующим (6.9) многочленом. Аппроксимирующий многочлен не проходит через все узловые точки таблицы. Поэтому его степень m не зависит от числа узловых точек. При этом всегда m < n. Степень m может меняться в пределах . Если m=1, то мы аппроксимируем табличную функцию прямой линией. Такая задача называется линейной регрессией. Если m=2, то мы аппроксимируем табличную функцию квадратичной параболой. Такая задача называется квадратичной аппроксимацией. Если m=3, то мы аппроксимируем табличную функцию кубической параболой. Такая задача называется кубической аппроксимацией. Уточним метод наименьших квадратов: для табличной функции, полученной в результате эксперимента, построить аппроксимирующий многочлен (6.9) степени m, для которого сумма квадратов уклонений по всем узловым точкам минимальна, т.е. (6.10) Изменим вид многочлена Pm. Поставим на последнее место слагаемые, содержащие xm. На предпоследнее - слагаемые, содержащие xm-1 и т.д. В результате получим: (6.11) или При этом изменим индексы коэффициентов многочлена. Тогда условие (6.8) будет иметь вид: где xi и yi - координаты узловых точек таблицы, -неизвестные коэффициенты многочлена (6.11). aj, Необходимым условием существования минимума функции S является равенство нулю ее частных производных по каждой aj. 49 В результате получили систему линейных уравнений. Раскрывая скобки и перенося свободные члены в правой части уравнений, получим в нормальной форме систему линейных уравнений: (6.12) где aj - неизвестные системы линейных уравнений (6.12), - коэффициенты системы линейных уравнений (6.12), - свободные члены системы линейных уравнений (6.12), Порядок системы равен m+1. Программирование метода наименьших квадратов (МНК) Изменим индексацию в системе (6.12). В результате получим: (6.13) где - неизвестные системы линейных уравнений (6.13), 50 - коэффициенты системы линейных уравнений (6.13), - свободные члены системы линейных уравнений (11.13), (xi, yi) - координаты узловых точек табличной функции, N - количество узловых точек, m - степень аппроксимирующего многочлена вида: , (6.14) Алгоритм задачи: 1. Строим систему линейных уравнений (6.13). Определяем коэффициенты ck,j и свободные члены dk. Т.к. система (6.13) симметрична относительно главной диагонали, то достаточно определить только наддиагональные элементы системы. 2. Решаем систему (6.13) методом Гаусса. Находим коэффициенты aj многочлена (6.14). 3. Строим аппроксимирующий многочлен (6.14) и определяем его значение в каждой узловой точке Pi = Pm(xi). 4. Находим уклонение каждой узловой точки . 5. Находим сумму квадратов уклонений по всем узловым точкам . 6. Находим остаточную дисперсию . Для построения аппроксимирующего многочлена (6.11) и вычисления его значения в каждой узловой точке используем рациональную форму многочлена: (6.15) Тогда для вычисления значения многочлена (6.15) удобно пользоваться схемой Горнера. Рекуррентная формула по схеме Горнера имеет вид: Укрупненная схема алгоритма МНК представлена на рис.6.7. 51 Рис. 6.7. Укрупненная схема алгоритма аппроксимации методом наименьших квадратов 52 Лекция 7. Компьютерное моделирование и решение дифференциальных уравнений В лекции рассматриваются методы моделирования систем, в которых входные переменные являются функциями от времени или каких-либо других параметров. Динамические системы - это системы, в которых входные переменные являются функциями от времени или каких-либо других параметров. Описываются эти системы дифференциальными или интегральными уравнениями. Например, большая часть законов механики, электротехники, теории упругости, теории управления и т.д. описываются с помощью дифференциальных уравнений. На практике динамические системы встречаются очень часто. Моделирование систем, связанных с движением тел, с расчетом потоков энергии, с расчетом потоков материальных ресурсов, с расчетом оборотов денежных средств и т.д. в конечном счете, сводится к построению и решению дифференциальных уравнений (как правило, II-го порядка). Прямолинейное движение тела, движущегося под действием переменной силы ,где S=S(t), описывается дифференциальным уравнением второго порядка в форме уравнения Ньютона: где m - масса тела, S - перемещение тела, -линейная скорость, -линейное ускорение. При этом задаваемые начальные условия имеют четкий физический смысл. Это - начальное положение тела и его начальная скорость. Вращательное движение тела под действием крутящего момента , где , описывается аналогично Где Iр - полярный момент инерции тела, 53 -угол поворота, - угловая скорость, - угловое ускорение. При построении математических моделей систем, машин, механизмов с учетом колебаний, возникающих в них, также необходимо построить и решить дифференциальное уравнение, т.к. все виды колебаний (свободные гармонические, вынужденные) также описываются дифференциальными уравнениями. На практике лишь небольшое число дифференциальных уравнений допускает интегрирование в квадратурах. Еще реже удается получить решение в элементарных функциях. Поэтому большое распространение при решении математических моделей с помощью ЭВМ получили численные методы решения дифференциальных уравнений. Нахождение определенного интеграла в процессе моделирования объектов процессов или систем может применяться в следующих задачах: 1. Определение пути при переменной скорости: 2. Нахождение скорости при переменном ускорении: 3. Определение моментов инерции тел: 4. Нахождение работы переменной силы: 5. При решении дифференциальных уравнений. Итак, дана функция y=f(x). Найти интеграл этой функции на участке [a,b], т.е. найти Если подынтегральная функция f(x) задана в аналитическом виде; если функция f(x) непрерывна на отрезке [a,b] ; если известна ее первообразная, т.е. 54 то интеграл может быть вычислен по формуле Ньютона-Лейбница как приращение первообразной на участке [a,b], т.е. Но на практике формула Ньютона-Лейбница для вычисления интеграла используется редко. Численные методы интегрирования применяются в следующих случаях: 1. подынтегральная функция f(x) задана таблично на участке [a,b] ; 2. подынтегральная функция f(x) задана аналитически, но ее первообразная не выражается через элементарные функции; 3. подынтегральная функция f(x) задана аналитически, имеет первообразную, но ее определение слишком сложно. В численных методах интегрирования не используется нахождение первообразной. Основу алгоритма численных методов интегрирования составляет геометрический смысл определенного интеграла. Интеграл численно равен площади S криволинейной трапеции, расположенной под подынтегральной кривой f(x) на участке [a,b] (рис.7.1). Рис. 7.1. Геометрический смысл определенного интеграла Суть всех численных методов интегрирования состоит в приближенном вычислении указанной площади. Поэтому все численные методы являются приближенными. При вычислении интеграла подынтегральная функция f(x) аппроксимируется интерполяционным многочленом. На практике чтобы не иметь дело с многочленами высоких степеней, весь участок [a,b] делят на части и интерполяционные многочлены строят для каждой части деления. Порядок вычисления интеграла численными методами следующий (рис.7.2): 1. Весь участок [a,b] делим на n равных частей с шагом h=(b-a)/n. 55 2. В каждой части деления подынтегральную функцию f(x) аппроксимируем интерполяционным многочленом. Степень многочлена n = 0,1,2: 3. Для каждой части деления определяем площадь частичной криволинейной трапеции. 4. Суммируем эти площади. Приближенное значение интеграла I равно сумме площадей частичных трапеций Рис. 7.2. Вычисление определенного интеграла Нахождение приближенного значения интеграла называется квадратурой, а формулы для приближенного вычисления интеграла - квадратурными формулами или квадратурными суммами. Разность R между точным значением интеграла и приближенным значением называется остаточным членом или погрешностью квадратурной формулы, т.е. Если в каждой из частей деления интервала [a,b] подынтегральная функция аппроксимируется многочленом нулевой степени, т.е. прямой, параллельной оси OX, то квадратурная формула называется формулой прямоугольников, а метод методом прямоугольников. Если в каждой из частей деления интервала [a,b] подынтегральная функция аппроксимируется многочленом первой степени, т.е. прямой, соединяющей две соседние узловые точки, то квадратурная формула называется формулой трапеций, а метод - методом трапеций. Если в каждой из частей деления интервала [a,b] подынтегральная функция аппроксимируется многочленом второй степени, то квадратурная формула называется формулой Симпсона, а метод - методом Симпсона. Метод прямоугольников 56 Словесный алгоритм метода прямоугольников: 1. Весь участок [a,b] делим на n равных частей с шагом h=(b-a)/n. 2. Определяем значение yi подынтегральной функции f(x) в каждой части деления, т.е. 3. В каждой части деления подынтегральную функцию f(x) аппроксимируем интерполяционным многочленом степени n = 0, т.е. прямой, параллельной оси OX. В результате вся подынтегральная функция на участке [a,b] аппроксимируется ломаной линией. 4. Для каждой части деления определяем площадь Si частичного прямоугольника. 5. Суммируем эти площади. Приближенное значение интеграла I равно сумме площадей частичных прямоугольников. Если высота каждого частичного прямоугольника равна значению подынтегральной функции в левых концах каждого шага, то метод называется методом левых прямоугольников (рис.7.3). Тогда квадратурная формула имеет вид Рис. 7.3. Метод левых прямоугольников Если высота каждого частичного прямоугольника равна значению подынтегральной функции в правых концах каждого шага, то метод называется методом правых прямоугольников (рис.7.4). Тогда квадратурная формула имеет вид 57 Рис. 7.4. Метод правых прямоугольников Точность каждого метода прямоугольников имеет порядок h. Алгоритм вычисления интеграла построим в виде итерационного процесса поиска с автоматическим выбором шага. На каждом шаге будем уменьшать шаг в два раза, то есть увеличивать число шагов n в два раза. Выход из процесса поиска организуем по точности вычисления интеграла. Начальное число шагов n=2.Схема алгоритма методов прямоугольников представлена на рис.7.5. Рис. 7.5. Схема алгоритма метода прямоугольников (с автоматическим выбором шага) Условные обозначения: 58 a,b - концы интервала, - заданная точность, с=0 - метод левых прямоугольников, с=1 - метод правых прямоугольников, S1 - значение интеграла на предыдущем шаге, S - значение интеграла на текущем шаге. Метод трапеций Словесный алгоритм метода трапеций: 1. Интервал [a,b] делим на n равных частей с шагом h=(b-a)/n. 2. Вычисляем значение подынтегральной функции в каждой узловой точке 3. На каждом шаге подынтегральную функцию f(x) аппроксимируем прямой, соединяющей две соседние узловые точки. В результате вся подынтегральная функция на участке [a,b] заменяется ломаной линией проходящей через все узловые точки. 4. Вычисляем площадь каждой частичной трапеции. 5. Приближенное значение интеграла равно сумме площадей частичных трапеций, т.е. Найдем площади Si частичных трапеций: Приближенное значение интеграла равно Точность метода трапеций имеет порядок h2. Схема алгоритма метода трапеций представлена на рис.12.6. 59 Рис. 7.6. Схема алгоритма метода трапеций (с автоматическим выбором шага) Метод Симпсона В методе Симпсона в каждой части деления подынтегральная функция аппроксимируется квадратичной параболой a0x2+a1x+a2. В результате вся кривая подынтегральной функции на участке [a,b] заменяется кусочно-непрерывной линией, состоящей из отрезков квадратичных парабол. Приближенное значение интеграла I равно сумме площадей под квадратичными параболами. Т.к. для построения квадратичной параболы необходимо иметь три точки, то каждая часть деления в методе Симпсона включает два шага, т.е. Lk=2h. В результате количество частей деления N2=n/2. Тогда n в методе Симпсона всегда четное число. Определим площадь S1 на участке [x0, x2] (рис.7.2). 60 Исходя из геометрического смысла определенного интеграла, площадь S1 равна определенному интегралу от квадратичной параболы на участке [x0, x2]: Неизвестные коэффициенты квадратичной параболы а0 , а1, а2 определяем из условия прохождения параболой через три узловых точки с координатами (x0y0), (x1y1), (x2y2). На основании этого условия строим систему линейных уравнений: Решая эту систему, найдем коэффициенты параболы. В результате имеем: .. Для участка [x2, x4]: ::::::::::::::::::: .. Для участка [xi-1, xi+1]: ., где . Суммируя все площади S1 под квадратичными параболами, получим квадратурную формулу по методу Симпсона: где N2 - количество частей деления. Точность метода Симпсона имеет порядок (h3/h4). Схема алгоритма метода Симпсона представлена на рис.7.7. 61 Рис. 7.7. Схема алгоритма Симпсона (с автоматическим выбором шага) Численные методы решения дифференциальных уравнений первого порядка Общий вид дифференциального уравнения (7.1) Нормальная форма дифференциального уравнения (7.2) где y=y(x) -неизвестная функция, подлежащая определению, f(x,y) - правая часть дифференциального уравнения в нормальной форме, равная первой производной функции y(x). В функцию f(x,y) помимо аргумента x входит и сама неизвестная функция y(x). 62 Пример: - общий вид дифференциального уравнения первого порядка, - нормальная форма этого же уравнения. Если неизвестная функция у зависит от одного аргумента x, то дифференциальное уравнение вида называется обыкновенным дифференциальным уравнением. Если функция у зависит от нескольких аргументов, то такое дифференциальное уравнение называется дифференциальным уравнением в частных производных. Общим решением обыкновенного дифференциального уравнения является семейство функций у=у(х,с) (рис 7.8): Рис. 7.8. При решении прикладных задач ищут частные решения дифференциальных уравнений. Выделение частного решения из семейства общих решений осуществляется с помощью задания начальных условий: (7.3) т.е. начальной точки с координатами (х0, у0). Нахождение частного решения дифференциального уравнения 12.4) удовлетворяющего начальному условию (7.3) называется задачей Коши. В численных методах задача Коши ставится следующим образом: найти табличную функцию которая 63 удовлетворяет заданному дифференциальному уравнению (7.2) и начальному условию (7.3) на отрезке [a,b] с шагом h, то есть найти таблицу i x y x0 y0 1 x1 y1 2 x2 y2 3 x3 y3 .. ... ... n xn yn . Здесь h - шаг интегрирования дифференциального уравнения, a=x0 - начало участка интегрирования уравнения, b=xn - конец участка, n=(b-a)/h - число шагов интегрирования уравнения. На графике (рис 7.9) решение задачи Коши численными методами представляется в виде совокупности узловых точек с координатами (xi ,yi), . Рис. 7.9. Методы Рунге - Кутта Наиболее эффективными и часто встречаемыми методами решениями задачи Коши являются методы Рунге - Кутта. Они основаны на аппроксимации искомой функции у(х) в пределах каждого шага многочленом, который получен при помощи разложения функции у(х) в окрестности шага h каждой i-ой точки в ряд Тейлора: (7.4) Усекая ряд Тейлора в различных точках и отбрасывая правые члены ряда, Рунге и Кутта получали различные методы для определения значений функции у(х) 64 в каждой узловой точке. Точность каждого метода определяется отброшенными членами ряда. Метод Рунге - Кутта 1-го порядка (метод Эйлера) Отбросим в (7.4) члены ряда, содержащие h2, h3, h4:. Тогда . Так как Получим формулу Эйлера: (7.5) Так как точность методов Рунге-Кутта определяется отброшенными членами ряда (7.4), то точность метода Эйлера на каждом шаге составляет . Алгоритм метода Эйлера можно построить в виде двух программных модулей: основной программы и подпрограммы ELER, реализующей метод. Рис. 7.10. Схема алгоритма метода Эйлера Здесь (x,y) -при вводе начальная точка, далее текущие значения табличной функции, h -шаг интегрирования дифференциального уравнения, b -конец интервала интегрирования. Рассмотрим геометрический смысл метода Эйлера. 65 Формула Эйлера имеет вид: где . Тогда формула Эйлера принимает вид: где - тангенс угла наклона касательной к искомой функции у(x) в начальной точке каждого шага. Рис. 6.11. Геометрический смысл метода Эйлера В результате в методе Эйлера на графике (рис 7.10) вся искомая функция y(x) на участке [a,b] аппроксимируется ломаной линией, каждый отрезок которой на шаге h линейно аппроксимирует искомую функцию. Поэтому метод Эйлера получил еще название метода ломаных. В методе Эйлера наклон касательной в пределах каждого шага считается постоянным и равным значению производной в начальной точке шага xi. В действительности производная, а, значит, и тангенс угла наклона касательной к кривой y(x) в пределах каждого шага меняется. Поэтому в точке xi+h наклон касательной не должен быть равен наклону в точке xi. Следовательно, на каждом шаге вносится погрешность. Первый отрезок ломаной действительно касается искомой интегральной кривой y(x) в точке (x0,y0). На последовательных же шагах касательные проводятся 66 из точек (xi,yi), подсчитанных с погрешностью. В результате с каждым шагом ошибки накапливаются. Основной недостаток метода Эйлера - систематическое накопление ошибок. Поэтому метод Эйлера рекомендуется применять для решения дифференциальных уравнений при малых значениях шага интегрирования h. Метод Рунге - Кутта 2-го порядка (модифицированный метод Эйлера) Отбросим в (7.4) члены ряда, содержащие h3, h4, h5:. Тогда (7.6) Чтобы сохранить член ряда, содержащий h2, надо определить вторую производную y"(xi).Ее можно аппроксимировать разделенной разностью 2-го порядка Подставляя это выражение в (7.6), получим Окончательно, модифицированная или уточненная формула Эйлера имеет вид: (7.7) Как видно, для определения функции y(x) в точке i+1 необходимо знать значение правой части дифференциального уравнения f(xi+1, yi+1) в этой точке, для определения которой необходимо знать предварительное значение yi+1. Для определения предварительного значения yi+1 воспользуемся формулой Эйлера. Тогда все вычисления на каждом шаге по модифицированной или уточненной формуле Эйлера будем выполнять в два этапа: На первом этапе вычисляем предварительное значение Эйлера по формуле На втором этапе уточняем значение y=i+1 по модифицированной или уточненной формуле Эйлера 67 Точность метода определяется отброшенными членами ряда Тейлора (12.4), т.е. точность уточненного или модифицированного метода Эйлера на каждом шаге . Рассмотрим геометрический смысл модифицированного метода Эйлера. Так как то модифицированную формулу Эйлера можно представить в виде: где - тангенс угла наклона касательной к искомой функции у(х) в начальной точке каждого шага, - тангенс угла наклона касательной к искомой функции у(х) в конечной точке каждого шага. Рис. 7.12. Геометрический смысл модифицированного метода Эйлера Здесь P1 - накопленная ошибка в (i+1)й точке по методу Эйлера, P2 - накопленная ошибка в (i+1)й точке по модифицированному методу Эйлера. 68 Как видно из рис.7.11, в первой половине каждого шага, то есть на участке [xi, xi+h/2], искомая функция y(x) аппроксимируется прямой, которая выходит из . точки (xi, yi) под углом, тангенс которого Во второй половине этого же шага, т.е. на участке [xi + h/2,xi + h], искомая функция y(x) аппроксимируется прямой, которая выходит из точки с координатами под углом, тангенс которого В результате в модифицированном методе Эйлера функция у(х) на каждом шаге аппроксимируется не одной прямой, а двумя. Алгоритм модифицированного метода Эйлера можно построить в виде двух программных модулей: основной программы и подпрограммы МELER, реализующей метод (рис. 7.13). Рис. 7.13. Схема алгоритма модифицированного метода Эйлера Здесь (x,y) -при вводе начальная точка, далее текущие значения табличной функции, 69 h -шаг интегрирования дифференциального уравнения, b -конец интервала интегрирования. Метод Рунге - Кутта 4-го порядка Самое большое распространение из всех численных методов решения дифференциальных уравнений с помощью ЭВМ получил метод Рунге-Кутта 4-го порядка. В литературе он известен как метод Рунге-Кутта. В этом методе на каждом шаге интегрирования дифференциальных уравнений искомая функция y(x) аппроксимируется рядом Тейлора (7.4), содержащим члены ряда с h4: В результате ошибка на каждом шаге имеет порядок h5. Для сохранения членов ряда, содержащих h2,h3,h4 необходимо определить вторую y", третью y"' и четвертую y(4) производные функции y(x). Эти производные аппроксимируем разделенными разностями второго, третьего и четвертого порядков соответственно. В результате для получения значения функции yi+1 по методу Рунге-Кутта выполняется следующая последовательность вычислительных операций: Вывод формулы не приведен. Предоставляется возможность вывод формул выполнить самостоятельно. Алгоритм метода Рунге-Кутта (4-го порядка) можно построить в виде двух программных модулей: основной программы и подпрограммы Rk4, реализующей метод (рис 7.14). 70 Рис. 7.14. Схема алгоритма метода Рунге-Кутта 4-го порядка. Здесь (x,y) -при вводе начальная точка, далее текущие значения табличной функции, h -шаг интегрирования дифференциального уравнения, b -конец интервала интегрирования. Решение дифференциальных уравнений высоких порядков Методы Рунге-Кутта можно использовать не только для решения дифференциальных уравнений первого порядка но и для решения дифференциальных уравнений более высоких порядков Любое дифференциальное уравнение m-го порядка (7.8) можно свести к системе, состоящей из m уравнений первого порядка при помощи замен. Заменим: 71 В результате дифференциальное уравнение m -го порядка (7.8) сводится к системе, состоящей из m дифференциальных уравнений первого порядка: (7.9) Решением системы (7.2), а значит и дифференциального уравнения m -го порядка (7.1) является m табличных функций Решение дифференциальных уравнений второго порядка В задачах моделирования динамических систем наиболее часто приходится решать дифференциальные уравнения второго порядка. Общий вид дифференциальных уравнений второго порядка: (7.10) Нормальная форма дифференциальных уравнений второго порядка: (7.11) Пример Уравнение в общем виде Его нормальная форма (7.12) Дифференциальное уравнение второго порядка (7.11) можно свести к системе, состоящей из двух дифференциальных уравнений первого порядка при помощи замен. Заменим y1=y', Тогда y'1=y". В результате уравнение (7.11) сводится к системе, состоящей из двух дифференциальных уравнений первого порядка: 72 Для примера (7.12) эта система имеет вид: (7.13) Решением этой системы являются две функции y(x) и y1(x), где Сформулируем задачу Коши для системы, состоящей из двух дифференциальных уравнений второго порядка. Дана система (7.14) Даны два начальных условия: Необходимо проинтегрировать систему на участке [a, b] с шагом h. В численных методах задача Коши для системы (7.14) ставится следующим образом: Найти табличные функции и , т.е. найти таблицу i x y y1 x0 y0 (y1)0 1 x1 y1 (y1)1 2 x2 y2 (y1)2 3 x3 y3 (y1)3 ... ... ... ... n xn yn (y1)n Здесь h - шаг интегрирования дифференциального уравнения, a=x0 - начало участка интегрирования уравнения, b=xn - конец участка, n=(b-a)/h - число шагов интегрирования уравнения. На графике решением задачи Коши для системы, состоящей из двух дифференциальных уравнений первого порядка, является совокупность узловых точек (рис. 7.15). 73 При этом на каждом шаге, т.е. для каждого значения xi решением являются две узловые точки с координатами (xi, yi), (xi, (y1)i). Рис. 7.15. Для решения системы дифференциальных уравнений используем те же методы, что и для решения одного дифференциального уравнения первого порядка. При этом необходимо соблюдать условие: на каждом шаге интегрирования, т. е. в точках с координатами х1 , х2 , х3 , : , хn все уравнения системы надо решать параллельно. Для вычисления правых частей уравнений системы (7.14) необходимо сформировать подпрограмму PRAV. Вернемся к примеру (7.13). Здесь на каждом шаге в подпрограмме PRAV будем вычислять правые части каждого уравнения системы: Схема алгоритма решения системы (7.13) представлена на рис 7.16. 74 Рис. 7.16. Схема алгоритма решения системы (7.6) Здесь h - шаг интегрирования дифференциального уравнения, b - конец участка, n - число шагов интегрирования уравнения, x, y, y1 - при вводе начальные значения, далее - текущие значения табличной функции. Решение дифференциальных уравнений m-го порядка методом РунгеКутта (4-го порядка) Как уже было сказано, любое дифференциальное уравнение m-го порядка (7.15) сводится к системе, состоящей из m дифференциальных уравнений 1-го порядка 75 (7.6) Численным решением системы (7.9), а значит и дифференциального уравнения m-го порядка (7.8) является m табличных функций т.е. функция y(x) и все ее производные, включая производную (m-1)-го порядка. При этом каждая из табличных функций определяется на промежутке [a, b] с шагом h и включает n узловых точек. Таким образом, численным решением уравнения (7.8) или системы (7.9) является матрица порядка , (табл. 7.8) Таблица 7.1. i x y y1=y', y1=y''1, : y_m-1=y^(m-1) x0 y0 (y1)0 (y2)0 : (ym-1)0 1 x1 y1 (y1)1 (y2)1 : (ym-1)1 2 x2 y2 (y1)2 (y2)2 : (ym-1)2 3 x3 y3 (y1)3 (y2)3 : (ym-1)3 : : : : : : : n xn yn (y1)n (y2)n : (ym-1)n где m - порядок дифференциального уравнения, равен количеству столбцов матрицы, n = (b-a)/h - количество шагов интегрирования, равно количеству строк матрицы. Каждый j -й столбец матрицы - это массив решений одной j -й табличной функции по всем n шагам интегрирования. Каждая i -ая строка матрицы - это массив решений m табличных функций на одном i -ом шаге интегрирования. На графике решением дифференциального уравнения m -го порядка (12.8) является совокупность узловых интегрирования, т.е. каждому значению xi, точек скоординатами 76 точек. При этом каждому шагу , соответствуют m узловых При построении алгоритма задачи будем как и ранее расчет вести по шагам интегрирования, т.е. в цикле по При этом, как и ранее, на каждом i -ом шаге цикла будем рассчитывать решение дифференциального уравнения и тут же его печатать. Тогда нет необходимости формировать матрицу решений, а можно ограничиться формированием массива решений (7.15), который соответствует одной i-й строке матрицы: (7.17) где (7.18) Тогда при построении системы дифференциальных уравнений изменится индексация. Рассмотрим пример 7.19). Дано дифференциальное уравнение второго порядка (7.19) С учетом обозначений (7.18) имеем: Тогда дифференциальное уравнение (7.19) сводится к системе: (7.20) Вычисление правых частей уравнений системы (7.20) будем выполнять в подпрограмме PRAV. При этом подпрограмма PRAV будет иметь вид Обращение к подпрограмме PRAV(m, x, Y, F), где m -порядок системы, x - значение x на i -м шаге интегрирования, Y=[y(1), y(2), y(3), : , y(m)] - входной массив длиной m, F=[f(1), f(2),:, f(m)] - результат, массив значений правых частей уравнений системы (12.10) длиной m. 77 Программу решения системы дифференциальных уравнений реализуем в виде 3-х программных модулей: 1. основной программы, в которой организуем циклический процесс по всем шагам интегрирования, ; 2. подпрограммы RGK, в которой на каждом i -м шаге реализуется метод Рунге Кутта (4-го порядка) для системы дифференциальных уравнений m -го порядка. Здесь на каждом шаге интегрирования вычисляется массив решений Y длиной m. Для этого внутри подпрограммы RGK организуем циклы по j, где ; 3. подпрограммы PRAV, обращение к которой осуществляется из подпрограммы RGK для вычисления массива F длиной m - значений правых частей уравнений системы. Схема алгоритма основной программы представлена на рис.7.17. Рис. 7.17. Схема алгоритма основной программы Здесь m -порядок системы, h -шаг интегрирования, n -количество шагов интегрирования, x -начальное и далее - текущее значение x, Y -массив длиной m, куда заносим начальные и далее - текущие значения решений системы на одном шаге интегрирования. В подпрограмме RGK для вычисления элементов массива Y, используем те же формулы, что и для решения одного дифференциального уравнения 1-го порядка методом Рунге-Кутта (4-го порядка), но с учетом поправки на массивы. Тогда 78 Здесь Y - массив решений длиной m, Y1 - рабочий массив длиной m, T - рабочая матрица порядка ( ). Для вычисления значений fj правых частей дифференциальных уравнений системы перед вычислением каждого элемента матрицы Т на каждом j -ом шаге необходимо выполнить обращение к подпрограмме PRAV. Схема алгоритма подпрограммы RGK представлена на рис.7.18 Рис. 7.18. Схема алгоритма подпрограммы RGK 79 Лекция 8. Метод конечных элементов как инженерный метод расчета. Краткая история развития МКЭ МКЭ является численным методом решения дифференциальных уравнений, встречающихся в физике и технике. Возникновение этого метода связанно с решением задач космических исследований (1950г.), и первые он был опубликован в работе М.Тернера, Р.Клужа, Г.Мартина и Л.Топпа (Turner M.J., Clouhg R. W., Martin H.C., Topp L.J. Stiffness and Deflection Analysis of Complex Structures // J. Aeronaut. Sci. – 1956. – №23. – P.805-824). Эта работа способствовала появлению других работ – был опубликован ряд статей с применениями метода конечных элементов к задачам строительной механики и механики сплошных сред. Важный вклад в теоретическую разработку метода сделал в 1965г. Р.Мелош (Melosh R.J. Basis for Derivation of Matrices for the Direct Stiffness method. // J. Am. Inst. For Aeronautics and Astronautics. –1965. - №1. – P.1631-1637), после чего показанная им связь МКЭ с процедурой минимизации функционала привела к широкому использованию МКЭ при решении задач в других областях техники. В первых работах с помощью метода решались задачи распространения тепла. Затем МКЭ был применен к задачам гидромеханики. Область применения существенно расширилась, когда О. Зенкевичем (Зенкевич О. Метод конечных элементов в технике. – М.: Мир. – 1975. – 541с) на основе глубокого анализа развития и апробации метода было показано, что уравнения, определяющие элементы в задачах строительной механики, распространения тепла, гидромеханики, могут быть легко получены с помощью таких вариантов метода взвешенных невязок, как метод Галеркина и метод наименьших квадратов. Установление этого факта сыграло важную роль в теоретическом обосновании МКЭ, так как позволило применять его при решении любых дифференциальных уравнений. МКЭ из численной процедуры решения задач строительной механики превратился в общий метод численного решения дифференциального уравнения или системы дифференциальных уравнений, в том числе и краевых задач теории упругости и теории пластичности. В СССР большой вклад в развитие МКЭ и его применение к прочностным расчетам в машиностроении внес уфимский ученый Р.Р.Мавлютов (Мавлютов Р.Р. Концентрация напряжений в элементах конструкций. – М.: Наука. – 1996. – 240с), которым показано, что МКЭ является одним из наиболее эффективных методов расчета. Он, в частности, позволяет с высокой точностью описать геометрию деталей сложной конфигурации, их напряженно-деформированное состояние в 80 зонах больших градиентов напряжений. С помощью МКЭ не представляет затруднений расчет конструкций из разнородных материалов, просто и точно учитываются реальные граничные условия, характеризующие контактные взаимодействия, адгезионные эффекты и т.п. Впоследствии МКЭ был развит для расчета процессов больших пластических деформаций Бурное развитие технологии расчетов МКЭ и прогресс вычислительной техники позволили применить МКЭ для моделирования столкновений автомобилей как с целью выявления недостатков конструкции, отрицательно влияющих на безопасность пассажиров, так и реконструкции обстоятельств ДТП. В середине 80-х годов прошлого века в США были построены и проанализированы первые полные модели автомобильных аварий, а их промышленное применение стало возможным с появлением в то время первых суперкомпьютеров. Так, например, уже в то время правительство Германии финансировало проект для исследования возможности численного моделирования автомобильных аварий на примере двух моделей автомобилей – Фольксваген-Поло и БМВ-300. Развитие применения численных методов в последующие десятилетия привело к тому, что МКЭ сегодня является инструментом, полностью интегрированным в процесс проектирования транспортного средства и элементов дороги, обеспечивающих безопасность. Сейчас конкурентно способное развитие отрасли невозможно без МКЭ-систем проектирования, которые середины 90-х годов используется всеми ведущими уже с автомобилестроительными компаниями. Так, например, фирма Меседес-Бенц для всех важных случаев ударного нагружения располагает детальными конечноэлементными аналогами более 30 моделей автомобилей с числом элементов более 200 тысяч каждая, и моделями манекенов водителя и пассажиров, которые непрерывно модифицируются, чтобы отслеживать соответствие требованиям стойкости при авариях (Du Bois P.A. Crashworthiness Engineering with LS-DYNA. – Livermore Software Technology Corporation, 2002). В любом случае МКЭ является единственным выбором, так как многочисленные требования и стандарты безопасности превышают возможности организации и анализа результатов натурных краш-тестов. Рассмотрим некоторые из множества примеров проверки расчетов МКЭ для ударов автомобилей. 81 Рис.8.1. Расчетная фактическая деформации автомобиля Додж Неон. Национальным центром анализа аварий (NCAC) университета им.Д.Вашингтона (США) на конечно-элементном аналоге из более чем 270 тысяч элементов был смоделирован краш-тест автомобиля Додж Неон 1996 года, произведенный сертифицированной лабораторией «Центр исследования транспорта» из штата Огайо по контракту с Департаментом транспорта США. Скорость фронтального удара автомобиля в жесткий недеформируемый неподвижный барьер была 56км/ч. На рис.1 показано сравнение расчетной и фактической деформации автомобиля. Рис.8.2. Расчетное (штриховая линия) и фактическое (сплошная линия) замедление (слева), скорость (в центре) и перемещение (справа) центра масс автомобиля Додж в зависимости от времени. Как видно из рис.8.1, совпадение расчетной и фактической формы деформированного автомобиля очень хорошее. На рис.8.2 приведено 82 сопоставление расчетных и фактических замедлений, скоростей и перемещений центра масс. Из рис. 8.2 видно, что расхождение расчетных и фактических параметров удара незначительное, и, например, расчетное и фактическое конечное положение центра масс автомобиля различаются не более чем на 50мм. Тем же NCAC был всесторонне исследован (Zaouk A., Bedewi N., Kan C., Marzougui D. Validation of non-linear finite element vehicle model using multiple impact date. – The George Washington University, NCAC) МКЭ-аналог пикапа Шевроле С-1500, показанный на рис. 8.3. Рис. 8.3. МКЭ-аналог пикапа Шевроле С-1500. Сначала аналог была испытан на фронтальный удар в плоский неподвижный жесткий барьер. Как и аналог Доджа, аналог пикапа Шевроле показал хорошее совпадение с результатами краш-теста. Затем было произведено испытание на скользящий удар пикапа в бетонное дорожное ограждение на скорости около 100км/ч. На рис. 8.4 показано сопоставление кадров видеосъемки с расчетными результатами на виде сверху. Видно, что результат расчета хорошо согласуется с фактическим движением пикапа. 83 Рис. 8.4. Сравнение фактического (вверху) и расчетного (внизу) движения пикапа при скользящем ударе. Для полноты на рис. 8.5 показано сопоставление кадров видеосъемки с расчетными результатами на виде спереди. Рис. 8.5. Сравнение расчетного (вверху) и фактического (внизу) движения пикапа при скользящем ударе. Здесь следует отметить, что скользящий удар является длительным и, как правило, корректно не воспроизводится иными, чем расчет МКЭ, методами реконструкции ДТП. 84 Рис. 8.6. Расчетная (жирная линия) и фактическая (тонкая линия) зависимость скорости пикапа от времени. На рис. 8.6 показана расчетная и фактическая зависимость скорости центра масс пикапа от времени. Видно, что в период скольжения пикапа по ограждению скорости совпадают с высокой точностью. После отделения пикапа от ограждения расчетная скорость отличается от фактической скорости на величину 1-2м/с=3.6-7.2км/ч, причем расчетная скорость выше фактической. Это обусловлено тем, что расчетная величина затрат кинетической энергии на деформацию меньше фактической, так как каким бы ни был подробным МКЭаналог автомобиля, учесть все элементы конструкции не представляется возможным. Последнее обстоятельство, – неполнота модели, – является существенным для использования МКЭ для судебной экспертизы ДТП. Вычисленная величина затрат энергии на деформацию всегда будет меньше фактической величины энергии, что приведет к установлению возможной наименьшей скорости движения ТС перед ударом или началом торможения. Так, например, при испытании пикапа на фронтальный удар в жесткий барьер фактическая величина кинетической энергии, затраченной деформацию, составила 232КДж, а расчетная величина – только 214КДж. на Ниже приведена таблица, показывающая распределение энергии, затраченной на деформацию, по элементам конструкции пикапа, как при фронтальном, так и скользящем ударе. Из таблицы следует, что неполнота 85 модели, т.е. учет деформации только основных силовых элементов конструкции кузова, позволяет установить не менее 80% фактически затраченной на деформацию кинетической энергии. Так как скорость (или изменение скорости) автомобиля в момент столкновения пропорциональна квадратному корню из величины затрат энергии, отсюда следует, что исходя из результатов прочностного расчета всегда можно установить не менее 90% фактической скорости (изменения скорости) автомобиля в момент столкновения, а при учете расхода энергии на торможение погрешность определения скорости (изменения скорости) будет еще меньше. Часть конструкции Фронтальный удар Затрат Вся модель Колеса и шины Рама Бампер и его крепление Двигатель и его крепление Радиатор и рамкой Затрат Скользящий удар Затрат Затрат 214.0 93.2 26.1 23.0 21.8 100.0 43.55 12.20 10.75 147.75 29.28 26.78 8.40 16.11 15.19 100.00 19.82 18.13 5.69 10.90 10.28 15.2 10.19 4.76 3.22 Рис. 8.7. МКЭ-аналоги автомобиля Форд и его рамы. Данные, приведенные в таблице, подтверждаются и иными исследованиями. Так, например, в 1998 году на международной конференции по реконструкции аварий был представлен доклад (Kirkpatrick S.W., Simons J.W., Antoun T.H. Development and validation of high fidelity vehicle crash simulation 86 models. // IJCrash’98 – International Crashworthiness Conference), авторы которого провели сравнение фактической и расчетной силы фронтального удара автомобиля Форд и его рамы в жесткий плоский барьер на скорости 56 км/ч. При этом расчетная величина деформирующей силы для полного МКЭ-аналога хорошо совпадала с фактической, установленной в краш-тесте. На рис. 8.7 показаны МКЭмодели автомобиля и его рамы. На рис. 8.8 показано сравнение зависимостей силы фронтального удара всего автомобиля и его рамы от времени. Рис.8. Зависимость силы удара от времени для полной модели Форда (сплошная линия) и его рамы (пунктирная линия). Из рис. 8.8 видно, что сила сопротивления рамы деформированию меньше, что естественно, силы сопротивления всей конструкции автомобиля. Однако пиковые значения сил сопротивления различаются незначительно, а для небольших деформаций графики практически совпадают. Рис. 8.9. МКЭ-аналог и фотография передней левой части Мерседеса. 87 На рис. 8.9 показаны МКЭ-аналог силовых элементов передней части и фотография передней левой части Мерседеса при снятом бампере. На рис. 8.10 показаны расчетная зависимость силы сопротивления лонжерона от величины деформации и фактическая зависимость этой силы от времени, известная из краштеста, произведенного департаментом транспорта США. Видно, что результаты хорошо совпадают. Так как силы сопротивления силовых элементов конструкций ТС существенно превышают их массу и силы сцепления шин с дорогой, из последних двух примеров исследований следует, что результаты прочностных расчетов части конструкции ТС могут применяться в судебной экспертизе для установления движения ТС в результате удара. Рис. 8.10. Расчетная зависимость силы удара от деформации (слева) и фактическая зависимость силы удара от времени (справа) для лонжерона Мерседеса. Известны исследования (Consolazio G.R., Chung J.H., Gurley K.R. Impact simulation and full scale crash testing of a low profile concrete work zone barrier. // Computers and Structures. – 2003. – V.81. – P.1359-1374), которые показывают, что относительно полные трехмерные МКЭ-модели более чем успешно конкурируют с импульсными моделями ДТП, известными, в частности, в России по таким программам моделирования ДТП, как CARAT или PC-Crash. На рис. 8.11 показан фрагмент МКЭ-аналога пикапа, содержащий колесо и подвеску. 88 Рис. 8.11. МКЭ-аналог колеса и подвески пикапа. Остальная часть МКЭ модели пикапа может быть и достаточно грубой, но должна хорошо отражать реальное распределение массы автомобиля. На рис. 8.12 показано сопоставление фактического и расчетного движения автомобиля при скользящем упругом контакте колесом с дорожным ограждением. Вывод из анализа результатов исследований, показанных на рис. 8.12, как и на рис. 8.4 с рис. 8.5, однозначный – МКЭ-аналоги позволяют с высокой точностью анализировать движение ТС в результате удара и могут успешно применяться при реконструкции обстоятельств ДТП в рамках судебной экспертизы. Прочностныерасчеты успешно применяются не только для установления обстоятельств ДТП, железнодорожных, авиационных аварий или аварий водных судов, но и для выяснения причин разрушения препятствий, когда разрушение препятствий имеет или может иметь 89 Рис. 8.12. Фактическое и расчетное движение автомобиля при упругом скользящем контакте с дорожным ограждением. катастрофические последствия. К таким препятствиям, в частности, относятся опоры мостов и остановочные павильоны, прочность которых должна обеспечить безопасность при возможных наездах ТС, или скорость ТС должна быть ограничена до величины, обеспечивающей устойчивость этих конструкций. На рис. 8.13 показан результат удара грузовика с полуприцепом в среднюю опору моста, проходящего через шоссе. От экспертов (El-Tawil S., Severino E., Fonseca P. Vehicle Collision with Bridge Piers. // Journal of Bridge Engineering. – 2005. – May/June. – P.353) требовалось установить причины разрушения моста. 90 На рис. 8.14 показана реконструкция удара грузовика в опору моста на скорости 110км/ч, произведенная расчетом МКЭ. Рис. 8.13. Авария 23.05.2003г. на шоссе I-80, штат Небраска, США. Результат удара грузовика с полуприцепом в опору моста. Рис. 8.14. Реконструкция удара грузовика в опору моста. 91 Рис. 8.15. Деформация опоры моста в результате удара На рис. 8.15 показана расчетная деформация опоры моста в течение удара, в результате которого произошел сдвиг верхней части опоры и обрушение перекрытий. Рис. 8.16. Расчетная деформация грузовика и опоры безопасности. На рис. 8.16 показаны расчетные деформации грузовика и опоры безопасности остановочного павильона, полученные в результате наезда грузовика на скорости 48км/ч (Lan S., Crawford J.E., Xin X. Development of Shallow Footing Anti-ram Bollard System Through Modeling and Testing. // 1-st International Conference on Analysis and Design). 92 Рис. 8.17. Результаты натурного эксперимента. На рис. 8.17 показаны результаты экспериментального теста – фотографии грузовика и опор перед наездом, фотография грузовика после наезда, фотография опоры после наезда и сопоставление расчетной и фактической зависимостей скоростей грузовика от времени. Сравнение рис. 8.16 и рис. 8.17 показывает хорошее совпадение расчетных фактических данных. Строгая научность, проработанность, точность МКЭ делает этот метод идеальным инструментом для инженерных расчетов. 93