0 МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ Белорусский национальный технический университет Кафедра «Строительная механика» ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ СТРОИТЕЛЬСТВА Лабораторный практикум Минск БНТУ 2016 1 МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ Белорусский национальный технический университет Кафедра «Строительная механика» ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ СТРОИТЕЛЬСТВА Лабораторный практикум для студентов специальности 1-70 02 01 «Промышленное и гражданское строительство» Рекомендовано учебно-методическим объединением высших учебных заведений по образованию в области строительства и архитектуры в качестве лабораторного практикума Минск БНТУ 2016 2 УДК 69:519.6 ББК 22.193:38.112 Ч67 Составители : А. В. Стрелюхин, Г. С. Богомолова, Е. Л. Сорокина Рецензенты : Р. И. Фурунжиев, А. С. Кравчук Численные методы решения задач строительства: лаборатор- ный практикум для студентов специальности 1-70 02 01 «Промыш- ленное и гражданское строительство» / сост. : А. В. Стрелюхин, Г. С. Богомолова, Е. Л. Сорокина. – Минск : БНТУ, 2016. – 116 с. ISBN 978-985-550-512-0. Настоящий практикум составлен на основе рабочей программы по курсу «Чис- ленные методы решения задач строительства». Материал состоит из введения, вось- ми лабораторных работ по численным методам решения задач линейной алгебры, дифференциальных уравнений с начальными и краевыми условиями, задач оптими- зации, основам работы в системе компьютерной алгебры «Mathematica», списка литературы и задач для самостоятельного решения. УДК 69:519.6 ББК 22.193:38.112 ISBN 978-985-550-512-0 © Белорусский национальный технический университет, 2016 Ч67 3 СОДЕРЖАНИЕ Введение .................................................................................................. 4 Лабораторная работа № 1. Введение в систему компьютерной алгебры Mathematica .................. 5 Лабораторная работа № 2. Решение систем линейных алгебраических уравнений методом Гаусса ..................................................................................... 20 Лабораторная работа № 3. Решение систем линейных алгебраических уравнений итерационными методами ................................................................... 31 Лабораторная работа № 4. Нахождение наибольшего по модулю собственного значения и соответствующего ему собственного вектора квадратной матрицы ............................................................................. 43 Лабораторная работа № 5. Решение задачи Коши методами Рунге-Кутта ................................... 57 Лабораторная работа № 6. Решение краевой задачи для линейных дифференциальных уравнений методом конечных разностей ........................................... 73 Лабораторная работа № 7. Численные методы оптимизации. Графический метод решения задач линейного программирования ............................................................. 88 Лабораторная работа № 8. Симплекс-метод решения задач линейного программирования ...... 99 Литература .......................................................................................... 112 Приложение......................................................................................... 113 4 ВВЕДЕНИЕ Лабораторный практикум предназначен для студентов дневной и заочной формы обучения специальности 1-70 02 01 «Промышлен- ное и гражданское строительство». Успешное освоение учебных дисциплин: «Строительная меха- ника», «Теоретическая механика», «Теория упругости» и т. д., со- ставляющих основу знаний современного инженера-строителя, – требует комплексной подготовки, в том числе и владения матема- тическим аппаратом. Курс «Численные методы решения задач строительства» носит междисциплинарный базовый характер, поз- воляет применять полученные знания для расчета заданий по дру- гим учебным дисциплинам и решать инженерные и экономические задачи строительства. Практикум состоит из восьми лабораторных работ, тематика ко- торых охватывает основные разделы курса «Численные методы решения задач строительства». Каждый раздел начинается с изло- жения теоретических основ методов (решение задач линейной алгебры, дифференциальных уравнений с начальными и краевыми условиями и оптимизации), далее рассматриваются типовые приме- ры с краткими указаниями, поясняющими алгоритм решения. Дополнительную информацию по разделам можно получить из ли- тературы, список которой приведен в конце практикума. В настоящее время одним из основных направлений развития программного обеспечения является разработка систем компьютер- ной алгебры, позволяющих переложить многие громоздкие и тру- доемкие вычисления на компьютер и сосредоточиться на проблеме постановки задачи и анализе результатов ее решения. Среди таких систем широкое применение получила система компьютерной алгебры «Mathematica». В связи с этим в практикуме уделяется осо- бое внимание изучению основ работы с системой «Mathematica» и применению полученных знаний для решения конкретных задач в соответствии с программой курса. Так же в практикуме приводятся задачи для самостоятельного решения. 5 ЛАБОРАТОРНАЯ РАБОТА № 1 ВВЕДЕНИЕ В СИСТЕМУ КОМПЬЮТЕРНОЙ АЛГЕБРЫ MATHEMATICA Цель работы: изучить основные возможности системы компью- терной алгебры Mathematica и уметь применять ее для инженерных расчетов. Состав работы: 1. Краткое описание системы компьютерной алгебры Mathe- matica; 2. Работа со списками, векторами и матрицами; 3. Работа с графикой; 4. Функция пользователя; 5. Выполнение индивидуальных заданий. 1.1. Краткое описание системы компьютерной алгебры Mathematica Решение многих инженерных задач требует проведения доста- точно большого количества расчетов. Кроме того, очень часто при- ходится обрабатывать полученные результаты, а также представ- лять их в наглядном виде. Выполнение таких операций вручную является достаточно трудоемкой работой. В связи с этим был разра- ботан целый ряд программных продуктов (так называемые системы компьютерной алгебры Computer Algebra System (CAS)) для обра- ботки числовых, символьных и графических данных, которые инте- грируют в себе современный интерфейс пользователя, возможности решения задач, как аналитически, так и численно, а также мощные средства визуализации полученных результатов с использованием различных типов графиков. Примером такого программного продукта является система ком- пьютерной алгебры Mathematica, разработанная фирмой Wolfram Research Inc. Первая версия системы появилась в 1988 г. В 2013 г. вышла девятая версия системы, а в 2014 г. – десятая версия. Возможности системы: аналитические преобразования (решение систем полиномиальных и тригонометрических уравнений и нера- венств, трансцендентных уравнений, нахождение пределов, инте- грирование и дифференцирование функций, решение дифференци- 6 альных уравнений и уравнений в частных производных и т. д.); численные расчеты (решение систем уравнений, нахождение пре- делов, интегрирование и дифференцирование, решение дифферен- циальных уравнений и уравнений в частных производных и т. д.); линейная алгебра (операции с матрицами, поиск собственных зна- чений и собственных векторов); графика и звук и многое другое. Запуск системы на компьютере под управление операционной системы Windows (при стандартной установке) осуществляется следующим способом Пуск – Программы – Wolfram Mathematica – Wolfram Mathematica 9 Интерфейс системы состоит из нескольких окон. На рисунке показан интерфейс программы с открытыми окнами главного меню системы, ввода информации для вычислений и окна Basic Math Assistant. Кратко рассмотрим содержание отдельных пунктов главного меню: File – предназначен для созда- ния нового окна документа, откры- тия существующего на внешнем носителе файла, закрытия окна до- кумента, сохранения полученных результатов в файл, предваритель- ного просмотра документа и его печати; Edit – содержит стандартные команды для работы с буфером обмена, поиском информации в тексте документа и настройке системы; Palettes – позволяет вызвать палитры (окна с набором пикто- грамм), с помощью которых можно получить быстрый и удобный доступ к вводу математических символов и функций; Window – содержит стандартные команды для работы с окнами; Help – содержит команды для доступа к справочной информа- ции, например к Documentation Center и Function Navi- gator, а также доступ к Web узлу компании Wolfram. 7 Работа с пунктами меню и их подпунктами осуществляется ана- логично другим приложениям под управлением операционной системы Windows. Ввод информации для расчетов основан на концепции рабочего поля или окна документа, называемого блокнотом (notebook). В нем могут находиться вычисления, текст, графика, звук и т. д., что в дальнейшем может использоваться для создания отчетов, пре- зентаций и книг. Блокноту соответствует файл текстового формата с расширением .nb. Основу блокнота составляют ячейки (cells). Ячейки указаны правыми скобками различных форм (по умолчанию скобка квад- ратная ]) и могут содержать в себе другие ячейки. Наиболее простой способ работы с системой Mathematica – интерактивный. В блокноте ввод данных (передача информации от «пользователя» к «вычислительному ядру» системы) произво- дится нажатием клавиш [Enter] [расширенной клавиатуры] или [Shift-Enter]. Клавиша [Enter] [основной клавиатуры] дает возможность вво- дить несколько команд. Например, если мы хотим вычислить выра- жение 5 + 7  8, необхо- димо набрать это выражение в окне ввода и запустить текущий блок на выполнение нажатием клавиш [Shift-Enter]. В результате выполнения данной операции получим результат 61. Признаком выполнения расчета является горизонтальная черта. В начале ввода информации на месте черты появится новый блок, отмеченный новой квадратной скобкой, в котором и будет отобра- жаться вводимая информация. При вычислениях исходное выражение присваивается объекту In[n], а результат вычисления – объекту Out[n], где n – номер проводимого вычисления. Значение n можно использовать в даль- нейших вычислениях. Объекты In[n] и Out[n] пользователем не вводятся за исключением случая обращения к конкретному объ- екту In[n] или Out[n]. 8 Например, прибавим 2 к полу- ченному в предыдущем примере результату. Точка с запятой (;) в конце записанного выражения пре- дотвращает вывод результата, но само вычисление выполняется. Комментарий в блокноте записывается с использованием (* … *) Символы %, %% и %n возвращают результат последней, предпо- следней и в строке n операции. Математические выражения в системе Mathematica записы- ваются с помощью операторов и функций. Ниже перечислены основные операторы для выполнения арифметических действий Операция Знак Операция Знак скобки ( ) Деление / факториал ! Умножение * или [пробел] двойной факториал !! Сложение + Возведение в степень ^ Вычитание – Намного упрощают работу по подготовке документов кнопочные палитры, предназначенные для ввода математических знаков. Если это ок- но не отображено, его можно вызвать из главного меню Palettes – Basic Math Assistant. Вместе с тем, для большинства операций можно использовать при вводе соответствующий набор кла- виш, например 9 Операция Комбинация клавиш Резуль- тат Операция Комбинация клавиш Резуль- тат Степень x[Ctrl+6]5 x5 Корень квадратный [Ctrl+2]x x Дробь x[Ctrl+/]5 5 x Радикал [Ctrl+2]x[Ctrl+5]5 5 x Система Mathematica имеет достаточно большой набор пред- определенных констант и функций, встроенных в ядро системы Mathematica, которые можно использовать при вычислениях. Основные константы Обозначение Описание Pi константа, хранящая значение числа π = 3,14159... E константа, хранящая значение числа е = 2,71828... I мнимая единица Infinity или ∞ знак бесконечность используется при некоторых вычислениях, например, при вычислении пределов, сумм и интегралов Degree константа для преобразования радиан в градусы, которая численно равно количеству радиан в одном градусе Основные элементарные функции Обозначение Описание Sqrt[x] квадратный корень Abs[x] абсолютное значение действительного или модуль комплексного числа x Mod[n,m] остаток от деления n на m Exp[x] показательная функция еx Log[x] натуральный логарифм ln(x) Log[a,x] логарифм по основанию a (logax) Sign[x] возвращает –1, 0 или 1, если аргумент х соответственно отрицательный, нулевой или положительный Тригонометрические функции Sin[x] синус sin(x), x в радианах Sin[x Degree] синус sin(x), x в градусах Cos[x] косинус cos(x), x в радианах Tan[x] тангенс tg(x), x в радианах Cot[x] котангенс ctg(x), x в радианах 10 Обратные тригонометрические функции ArcSin[x] арксинус arcsin(x) ArcCos[x] арккосинус arccos(x) ArcTan[x] арктангенс arctg(x) ArcCot[x] арккотангенс arcctg(x) Гиперболические функции Sinh[x] гиперболический синус Cosh[x] гиперболический косинус Tanh[x] гиперболический тангенс Coth[x] гиперболический котангенс Функции для работы с комплексными числами Arg[z] аргумент комплексного числа z Im[z] мнимая часть комплексного числа z Re[z] действительная часть комплексного числа z Conjugate[z] комплексно-сопряженное с z число Sign[z] отношение z/|z| комплексного числа z Полный перечень встроенных функций можно найти в спра- вочной системе программы (главное меню Help – Function Navigator). Справку по конкретной функции можно получить из окна доку- мента, установив на ней курсор и нажав клавишу F1, или непосред- ственно в ячейке, например ?Plot вызов справки по указанной функции ??Plot вызов расширенной справки по указанной функции ?Plo* вызов всех функций, содержащих первые буквы Plo в сво- ем имени ?*Plo* вызов всех функций, содержащих буквы Plo в своем имени При наборе выражений в системе Mathematica существуют опре- деленные правила: – дробная и вещественная часть числа отделяется точкой; – запятая разделяет аргументы функции; – заглавные и прописные буквы различаются; – имена встроенных функций всегда начинаются с заглавной буквы, например Cos. Если имя функции составное, т. е. состоит из нескольких корней, то каждая составная часть в имени функции начинается с заглавной буквы, например ArcCos. 11 – использование скобок различных типов жестко регламентиро- вано: круглые скобки ( ) используются для группировки выраже- ний и изменения порядка вычисления выражений, т. е. как и при обычных математических вычислениях. Использование других типов скобок, как это принято при обычных математических вычислениях, не допускается. квадратные скобки [ ] являются основным признаком функ- ций и используются при записи функций, внутри которых записы- ваются аргументы, разделенные запятыми, например Cos[x]. фигурные скобки { } используются при работе с массивами, списками и матрицами. Проиллюстрируем полу- ченную информацию на при- мере вычисления выражения 10sin(π/2)(1 – 5cos(π)). Так как изначально система Mathematica разрабатывалась для проведения аналитических расчетов, то не все выражения по умол- чанию вычисляются численно. Например: Для того, чтобы получить результат в численном виде, суще- ствуют два способа. Первый способ обычно используется для простых выражений и заключается в представлении числа как числа с плавающей запятой. Вторым способом является использование специальной функ- ции, которая сообщает системе, что данное выражение необходимо вычислить численно. 12 Эта функция называется N и имеет формат записи N[выражение, точность]. Если точность не указы- вается, то вычисления проводятся по умолчанию в пять знаков после точки. Результат использования различных способов показан на рисунке. При вычислении сложных выражений часто оказывается неудоб- но записывать его полностью. Тогда при его составлении можно воспользоваться переменными, которым присваивают результат промежуточного вычисления или некоторое значение. Имя пере- менной всегда начинается с буквы. Удобством переменных является то, что ими можно воспользоваться в любой момент после того как она была инициализирована и в любом месте вычислений. Следует помнить, что значения переменных, после того как они были про- считаны, хранятся только в течение текущего сеанса вычислений. При следующем запуске системы их значе- ния автоматически не восстанавливаются – требуется заново запустить соответствую- щие ячейки на расчет. Пример записи составного выражения показан на рисунке (в нем от вывода про- межуточных результатов мы отказались с помощью символа ;). Рассмотреть все возможности системы Mathematica не представ- ляется возможным, поэтому остановимся только на основных вычислениях. Дополнительные возможности системы будут рас- смотрены в последующих лабораторных работах. 1.2. Работа со списками, векторами и матрицами Списки являются не только эффективным средством работы с выражениями в процессе численных и символьных вычислений в системе Mathematica, но и одним из наиболее фундаментальных способов структурирования данных множественного типа в виде массивов, матриц и др. Основные сведения о матрицах и действиях над ними находятся в Приложении. 13 В системе Mathematica списки задаются с использованием фигурных скобок {}. Например: {1,2,3} список из 3 целых чисел {1.5,2.3,3.7} список из 3 вещественных чисел {a,b,c} список из 3 символьных данных {{a,b},{c,d}} список, эквивалентный матрице a b c d     С помощью списков можно задавать более привычные типы сложных данных, например, векторы – одномерные массивы дан- ных, матрицы – двумерные массивы и многомерные массивы. При этом каждая группа элементов многомерных списков выделя- ется своей парой фигурных скобок, а в целом список выделяется общими скобками. Например, список {{1,2,3},{4,5,6},{7,8,9}} представ- ляет собой матрицу 3 порядка 1 2 3 4 5 6 7 8 9       . Списки можно составлять напрямую, задавая объекты в соответ- ствии с описанным синтаксисом. Списки могут быть объектами присваивания значений переменным, например: Списки характеризуются размером, который представляет собой произведение числа элементов списков по каждому направлению. Число направлений называют размерностью. Например, одномер- ный список является вектором и характеризуется числом элементов по единственному направлению. При этом вектор может быть век- тором-строкой и вектором-столбцом. Двумерный список представ- ляет матрицу, имеющую m строк и n столбцов. Ее размер mn. Ес- ли m = n, матрица называется квадратной. Особенно часто для генерации списков с элементами целого и вещественного типа используется функция Table, создающая таблицу-список. V:={1,2,3,4,5} . 14 Table[expr,{i,imin,imax}] генерирует список значений expr при изменении i от значения imin до imax с шагом 1 Table[expr,{i,imin,imax,istep}] генерирует список значений expr при изменении i от значения imin до imax с шагом istep Например: Для вывода элементов списка используются следующие функции: – TableForm[список]выполняет вывод с элементами списка список, расположенными в массиве прямоугольных элементов; – MatrixForm[список]выводит список в форме матрицы. 15 Вектора в системе Mathematica трактуются как линейные, т. е. одноуровневые списки; матрицы – как двухуровневые. Другими словами можно сказать, что матрица является списком списков. Рассмотрим основные функции, предназначенные для работы с векторами и матрицами. Обозначение Описание Det вычисляет детерминант квадратной матрицы Inverse вычисляет обратную матрицу для невырожденных квадратных матриц Transpose вычисляет транспонированную матрицу Dot[s1,s2] или s1.s2 вычисляет скалярное произведение векторов, произве- дение вектора и матрицы, а также произведение матриц s1*s2 вычисляет поэлементное произведение списков Cross[s1,s2] для векторов вычисляет векторное произведение +, – сумма и разность векторов и матриц Norm[s1] вычисляет норму списка Использование отдельных функций при работе с матрицами рассмотрим на следующих примерах. Определитель матрицы Det Умножение матрицы на вектор и матрицы на скаляр 16 Примечание: 1. Для ввода элементов матрицы можно исполь- зовать соответствующую пиктограмму из палитры Basic Math Assista nt. Для добавления нуж- ного количества строк и столбцов используются кнопки Add Row и Add Column. 2. Для выполнения поэлементных действий с матрицами исполь- зуются двойные квадратные скобки: m1[[1,1]]+m2[[2,3]]. 1.3. Работа с графикой Система Mathematica имеет в своем распоряжении не очень большой набор функций для построения графиков, но, несмотря на это, она позволяет строить практически все типы графиков. Благодаря большому количеству опций и директив для графических функций система позволяет строить достаточно продвинутые и сложные графики с красивым оформлением. 1.3.1. Построение графиков функций, заданных аналитически Все функции для построения двухмерных графиков оканчи- ваются на Plot, а трехмерных – на Plot3D. Plot[f(x),{x,xmin,xmax}]. Построение графика функции f(x) в декартовых координатах для x в пределах от xmin до xmax. Plot3D[f(x,y),{x,xmin,xmax},{y,ymin,ymax}]. Построение трехмерного графика функции f(x, y) в декартовых координатах при изменении аргументов функции от минимального и до максимального значения. Использование функций для построения графиков рассмотрим на следующих примерах. 17 Функция одной переменной Две функции одной переменной на одном графике Функция двух переменных Две функции двух переменных на одном графике Кроме того, в системе есть отдельные функции для построения графика функции, заданной в параметрической форме (Paramet- ricPlot, ParametricPlot3D), контурного графика (ContourPlot), графика плотности (DensityPlot) и т. д. Более подробную информация о них, а также об опциях для различного отображения графиков, можно получить из справочной системы программы. 1.3.2. Построение графиков функций, заданных таблично Далеко не все задачи можно решить в аналитическом виде. В этом случае на помощь приходят численные методы. Результатом таких вычислений являются табличные данные, которые хотелось бы отобразить в виде графиков. Система Mathematica имеет несколько встроенных функций, аналогичных рассмотренным ранее функциям. 18 Основные функции для построения графиков по точкам следующие: – ListPlot[список]. Если в качестве параметра передает- ся обычный список, то его значения будут взяты в качестве коорди- нат по у при построении графика, а координаты по x будут заданы из списка 1, 2, 3, ... . Если же в качестве параметра будет передана матрица, состоящая из двух строк, то данные из первого столбца будут взяты в качестве x-координат, а второго столбца, соответ- ственно, в качестве y-координат для точек. График отображается в виде точек. – ListLinePlot[список]. Аналогична функции List- Plot, но график отображается в виде линий – ListPlot3D[матрица]. По данным, представленным в матрице с тремя столбцами, строится поверхность, где столбцы задают значения х, у и z-координат. График функции, заданный таблично График функции, заданный таблично Координаты x заданы из списка 1, 2, 3, … . Размер точки установлен опцией PlotStyle Координаты x и y заданы как эле- менты матрицы. Размер точки уста- новлен опцией PlotStyle. Опция Mesh→All отображает точки 1.4. Функция пользователя Хотя в системе Mathematica около тысячи встроенных функций, часто пользователю может потребоваться создание своей собст- 19 венной функции. Процесс создания функции пользователя рассмот- рим на следующем примере. Здесь мы определили собственную функцию 2 10( ) 5 20 xy x x    , присвоив ей имя Myfun. В ней x является формальным параметром и обязательно должен содержать символ подчеркивания при обра- щении к нему. 1.5. Выполнение индивидуальных заданий Индивидуальное задание этой лабораторной работы подразуме- вает самостоятельное выполнение в системе Mathematica приведен- ных выше примеров. 20 ЛАБОРАТОРНАЯ РАБОТА № 2. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ МЕТОДОМ ГАУССА Цель работы: изучить метод Гаусса решения систем алгебраи- ческих уравнений. Состав работы: 1. Изучение метода Гаусса (метод исключения неизвестных); 2. Изучение метода Гаусса с частичным выбором ведущего элемента; 3. Пример решения системы линейных алгебраических уравнений; 4. Решение систем линейный алгебраических уравнений в систе- ме Mathematica; 5. Выполнение индивидуальных заданий. 2.1. Метод Гаусса (метод исключения неизвестных) Пусть задана система линейных алгебраических уравнений (СЛАУ) n-го порядка 11 1 12 2 1 1 21 1 22 2 2 2 1 1 2 2 ... ; ... ; ... , n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b              (2.1) где n ≥ 2. В матричном виде исходную систему можно представить как A·X = B (основные сведения о матрицах находятся в Приложении). Требуется найти решение системы (2.1), т. е. найти такие значе- ния неизвестных x1, x2, … , xn, которые при подстановке во все уравнения исходной системы обращают их в тождества. Из линей- ной алгебры известно, что система имеет единственное решение, если ее определитель не равен нулю, т. е. система является невыро- жденной. Метод Гаусса состоит из двух основных этапов (ходов): 1) прямого хода; 2) обратного хода. 21 Цель прямого хода – преобразовать заданную систему к верхне- му треугольному виду. Это значит, что первое уравнение новой системы сохранит свой прежний вид. Из второго уравнения исклю- чается неизвестное x1, из третьего – x1 и x2 и т. д. Таким образом, последнее уравнение будет содержать только одно неизвестное xn. Цель обратного хода – непосредственно определить значения неизвестных, причем вначале определяется xn, затем xn-1, … , x1. Рассмотрим прямой ход. Чтобы исключить x1 из второго уравне- ния, найдем множитель 21 11/m a a (предполагаем, что коэффици- ент a11, который в методе Гаусса называется первым ведущим эле- ментом, отличен от нуля). Затем первое уравнение системы умно- жим на полученный множитель и почленно вычтем из второго уравнения. При этом его коэффициенты и свободный член получат новые значения (1) 21 21 11 21 1121 11 0;aa a a m a a a      (1) 22 1222 ;a a a m  ……… (1) 2 12 ;n nna a a m  (1) 2 12b b b m  . В результате мы получили, что во втором уравнении отсутствует x1. Чтобы исключить x1 из третьего и последующих уравнений, надо найти множители m = a31/a11, … , m = an1/a11, умножить на них все то же первое уравнение и почленно вычесть его из третьего и последующих уравнений. В результате всех этих действий заданная система примет вид 11 1 12 2 1 1 (1) (1) (1) 222 2 2 (1) (1) (1) 22 ... ; ... ; ... , n n nn nn n nn a x a x a x b a x a x b a x a x b            где первый и второй ведущие элементы специально выделены. 22 Переходим к исключению x2 из третьего и последующих уравне- ний, полагая, что второй ведущий элемент (1)22a не равен нулю. Находя множители (1) (1)32 22/m a a , … , (1) (1)2 22/nm a a , умножаем на них (поочередно) второе уравнение и почленно вычитаем из тре- тьего, … , n-го уравнений. Аналогично проводится исключение x3, … , xn–1. Все действия, направленные на исключение неизвестного xi из j-го уравнения, можно представить тремя формулами: ( 1) ( 1)/i iji iim a a   ; ( ) ( 1) ( 1)i i i jl jl ila a a m    (l = i +1, … , n); (2.2) ( ) ( 1) ( 1)i i i j j ib b b m    . Проведя вычисления по формулам (2.2) при изменении i от 1 до n–1, а j от i+1 до n, получим систему верхнего треугольного вида 11 1 12 2 1 1 (1) (1) (1) 222 2 2 ( 1) ( 1) ... ; ... ; , n n nn n n nn n n a x a x a x b a x a x b a x b           в которой ни один из ведущих элементов не должен быть равен нулю. Рассмотрим обратный ход. В последнее уравнение входит только одно неизвестное xn. Поэтому оно определится как ( 1) ( 1)/n nn n nnx b a   . (2.3) Предпоследнее уравнение преобразованной системы имеет вид ( 2) ( 2) ( 2) 11, 1 1, 1 n n n n nn n n n na x a x b        . 23 Подставляя в него найденное значение xn, получим  ( 2) ( 2) ( 2)1 1 1, 1, 1/n n nn nn n n n nx b a x a        . Продолжая этот процесс, достигнем первого уравнения, из кото- рого  1 1 12 2 1 11n nx b a x a x a    . Все действия, совершаемые на обратном ходе при определении xn-1, xn-2, … , x1, можно описать одной формулой ( 1) ( 1) ( 1) 1 / ni i i i i is s ii s i x b a x a          , (2.4) где i = n–1, … , 1. После окончания обратного хода найденные значения x1, x2, …, xn необходимо подставить во все уравнения исходной системы, чтобы убедится в том, что они обращаются в тождества и что x1, x2, …, xn образуют решение. Примечания. 1. Решая систему методом Гаусса, надо контролировать, чтобы какой-нибудь ведущий элемент не оказался равным нулю, так как деление на нуль невозможно, а именно в операции деления и участ- вуют ведущие элементы. Нулевой ведущий элемент может появится в двух случаях. Во-первых, когда система вырожденная и ее опре- делитель равен нулю. В этом случае решение вообще не существует либо имеется множество решений. Во-вторых, заданная система, будучи невырожденной, может иметь коэффициент a11 = 0. Для это- го используют, например, метод Гаусса с выбором ведущего элемента (раздел 2.2). 2. Описанный формулами (2.2) – (2.4) алгоритм можно обобщить на случай решения системы с k правыми частями. В этом случае свободные члены и неизвестные имеют два индекса и образуют матрицы (массивы) с размерами n  k. Третья из формул (2.2), а также формулы (2.3) и (2.4) теперь будут выглядеть так 24 ( ) ( 1) ( 1)i i i jl jl ilb b b m    ; ( 1) ( 1)/n nnl nnnlx b a   ; ( 1) ( 1) ( 1) 1 ( ) / ni i i il is sl iiil s i x b a x a        , где l = 1, … , k. 3. Метод Гаусса решения СЛАУ можно применить для вычисле- ния определителя матрицы А n-го порядка. В отличие от классиче- ской формулы, это выполняется очень просто, как произведение n ведущих элементов:    (1) ( 1)( 1)11 22 1 1 ... 1 n in nn ii i a a a a          , (2.5) где γ – число перестановок строк при выборе ведущего элемента в процессе преобразования исходной матрицы А к треугольно- му виду. 4. Метод Гаусса можно применить для поиска обратной матрицы. Исходя из определения обратной матрицы A·A-1 = E и используя для наглядности обозначение X = A-1, приходим к системе линей- ных алгебраических уравнений вида A·X = E, решение которой и будет составлять элементы обратной матрицы. 2.2. Метод Гаусса с частичным выбором ведущего элемента Пусть требуется решить систему трех уравнений 2 3 1 2 3 1 2 3 4,2 8,1 40,5, 2,2 4,1 6,3 5,6, 6,8 2,2 3,4 14,3. x x x x x x x x           Данная система является невырожденной (определитель равен – 397,07), следовательно имеет единственное решение. Но в данном случае метод Гаусса неприменим, т.к. первый ведущий эле- мент a11 = 0. 25 У метода Гаусса есть еще один существенный недостаток: веду- щий элемент может оказаться не точным нулем, а очень малой величиной. Тогда соответствующие множители получатся очень большими. Они внесут в преобразованные коэффициенты и свобод- ные члены большие погрешности, способные накапливаться в про- цессе дальнейшего решения. В таком случае полученные результа- ты могут оказаться неудовлетворительными. Выход из этих ситуаций заключается в выборе ведущего эле- мента. Самый простой вариант такого выбора состоит в том, что перед тем, как исключить x1 из второго и последующих уравнений, сравниваются между собой коэффициенты первого столбца a11, a21, … , an1, среди которых выбирается наибольший по модулю. Пусть таковым будет не a11, а некоторый другой – ak1. Тогда, поме- няв местами первое и k-е уравнения, на месте a11 получим величину, наибольшую из всех коэффициентов первого столбца. Вследствие этого соответствующие множители будут по модулю не большими, чем единица. В примере наибольший коэффициент первого столбца находится в третьем уравнении (k = 3). Поэтому меняем местами первое и третье уравнения: 1 2 3 1 2 3 2 3 6,8 2,2 3,4 14,3, 2,2 4,1 6,3 5,6, 4,2 8,1 40,5. x x x x x x x x           Определив множители m = 2,2/6,8 ≈ 0,324 и m = 0/6,8 = 0, преоб- разуем эту систему к виду 1 2 3 2 3 2 3 6,8 2,2 3,4 14,3 3,388 7,4 10,227 4,2 8,1 40,5. x x x x x x x          Перед исключением x2 из третьего и последующих уравнений выбор наибольшего по модулю ведется среди коэффициентов (1) (1) 22 2, ..., na a . В рассматриваемом примере наибольшим по модулю будет (1)32 4,2a   . Меняем местами второе и третье уравнения 26 1 2 3 2 3 2 3 6,8 2,2 3,4 14,3, 4,2 8,1 40,5, 3,388 7,4 10,227. x x x x x x x          Определив множитель m = –3,388/(–4,2) ≈ 0,807, приходим к системе 1 2 3 2 3 3 6,8 2,2 3,42 14,3, 4,2 8,1 40,5, 13,937 22,457. x x x x x x        Используя обратный ход, получаем x1 = 1,611; x2 = –6,536, x3 = 5,028. Определитель матрицы A определим с помощью выражения (2.5) det|A| = (–1)2·6,8·(–4,2)·13,937 ≈ –398,04. 2.3. Пример решения системы линейных алгебраических уравнений Решить СЛАУ третьего порядка методом Гаусса 1 2 3 1 2 3 1 2 3 25,2 4,7 1,4 46,8, 4,1 39,4 2,9 25,6, 1,3 3,1 16,8 34,3 x x x x x x x x x            Результаты представить с четырьмя значащими цифрами. Перед тем, как решать задачу, необходимо вычислить определи- тель системы: det|А| = 16884,793, т. е. матрица является невырож- денной и решение этой системы единственно. Выполним прямой ход. Исключаем x1 из второго уравнения. Для этого находим множи- тель 21 11/ 4,1 / 25,2 0,162698m a a   , первое уравнение системы умножаем на этот множитель и почленно вычитаем его из второго уравнения: 27 (1) 22 1222 39,4 ( 4,7) 0,162698 40,1647a a a m       , (1) 23 1323 2,9 1,4 0,162698 2,67222a a a m      , (1) 2 12 25,6 46,8 0,162698 33,2143b b b m        . Исключаем x1 из третьего уравнения аналогичным образом: 31 11/ 1,3 / 25,2 0,0515873m a a     , (1) 32 1232 3,1 ( 4,7) ( 0,0515873) 2,85754a a a m        , (1) 33 1333 16,8 1,4 ( 0,0515873) 16,8722a a a m       , (1) 3 13 34,3 46,8 ( 0,0515873) 36,7143b b b m       . В результате получили преобразованную систему: 1 2 3 2 3 2 3 25,2 4,7 1,4 46,8, 40,1647 2,67222 33,2143, 2,85753 16,8722 36,7143. x x x x x x x         Исключаем x2 из третьего уравнения: (1) (1) 32 22/ 2,85753 / 40,1647 0,0711453m a a   , (2) (1) (1) 33 33 23 16,8722 2,67222 0,0711453 16,6821a a a m      , (2) (1) (1) 3 3 2 36,7143 ( 33,2143) 0,0711453 39,0773b b b m       . Окончательно преобразованная система имеет вид 1 2 3 2 3 3 25,2 4,7 1,4 46,8, 40,1647 2,67222 33,2143, 16,6821 39,0773. x x x x x x        28 Теперь выполним обратный ход: (2) (2) 3 3 33/ 39,0773 /16,6821 2,34247,x b a   (1) (1) (1) 2 32 23 22( ) / ( 33,2143 2,67222 2,34247) / 40,1647 0,982801, x b a x a                1 1 12 2 13 3 11 46,8 4,7 0,982801 1,4 2,34247 25,2 1,54370. x b a x a x a            Проверка полученного решения. Подставим полученные значения x1, x2, x3 в левую часть исход- ной системы: 25,2 1,54370 4,7 ( 0,982801) 1,4 2,34247 46,7999,       4,1 1,54370 39,4 ( 0,982801) 2,9 2,34247 25,6000,        1,3 1,54370 3,1 ( 0,982801) 16,8 2,34247 34,3000        . Найденные числа с достаточно высокой точностью совпадают с правыми частями заданной системы. Это значит, что система решена правильно. Значения неизвестных округляем до четырех значащих цифр: x1 = 1,544, x2 = –0,9828, x3 = 2,343. Для того, чтобы вычислить определитель матрицы A по методу Гаусса, воспользуемся выражением (2.5): det|A| = (–1)0·25,2·40,1647·16,6821 ≈ 16 884,79. 2.4. Решение систем линейный алгебраических уравненийв системе Mathematica Для решения систем линейных уравнений в системе Mathematica используется функция LinearSolve[A,В], которая возвращает вектор X, являющийся решением матричного уравнения А·Х = В. Рассмотрим пример из п. 2.3 и воспользуемся системой Mathematica для решения. 29 Для ввода элементов мат- риц А и B использовалась пик- тограмма Matrix из палитры Basic Math Assistant (вкладка Matrix Commands). Нужное количества строк и столбцов устанавливается кнопками Add Row и Аdd Column 2.5. Выполнение индивидуальных заданий Составить систему линейных алгебраических уравнений, выбрав коэффициенты и свободные члены из табл. 2.1 согласно шифру (см. указания к работе). Решить эту систему методом исключения Гаусса. Сделать проверку полученных результатов. Решить задание, используя систему Mathematica, и сравнить полученные результаты. Таблица 2.1 Исходные данные к заданию Первая цифра шифра Коэффициенты системы уравнений Вторая цифра шифра Свободные члены bi Третья цифра шифра Вели- чина с ai1 ai2 ai3 1 2 3 4 5 6 7 8 0 4 + с 1,207 0,248 2,353 3,125 0,953 0,189 1,424 2,267 0 2,411 –0,315 2,000 0 0,954 1 3,742 4,025 1,562 6 + c 2,994 0,738 1,256 –0,984 4,092 1 –1,940 1,688 2,198 1 –1,025 30 Окончание табл. 2.1 1 2 3 4 5 6 7 8 2 3,116 –2,158 2,352 2,125 3,734 1,616 3 + c 3,450 3,528 2 1,480 2,914 –2,165 2 0,305 3 2,691 2 + c 3,546 2,472 3,924 4,017 –2,073 1,925 3,638 3 –1,270 –2,424 1,372 3 –0,167 4 3,886 3,004 –2,673 3,457 5 + c 2,392 1,968 3,785 4,101 4 1,801 –1,212 –2,624 4 0,425 5 5,116 2,158 3,195 2,125 4,168 –2,207 1,140 2 + c 4,218 5 –2,466 –1,215 1,272 5 –0,203 6 4,071 2,848 3 + c –3,797 2,442 4,504 3,376 1,955 3,973 6 2,837 –3,004 3,725 6 –0,632 7 5,113 2,957 3,306 2,588 4,625 4 + c –2,164 2,949 4,092 7 1,430 2,175 0,998 7 0,530 8 2,930 –3,468 2,504 –2,152 2,987 3,706 2,146 4,225 5 + c 8 –2,294 1,748 1,025 8 –0,758 9 6 + c 4,436 2,407 2,890 4,029 –3,913 2,476 –2,923 4,738 9 2,500 –1,692 1,573 9 0,873 Указания к работе. 1. Исходные данные к контрольной работе выбираются студен- том из таблицы в соответствии с его личным учебным шифром (номером зачетной книжки). Шифром считаются три последние цифры. Если номер зачетной книжки двузначный, то его следует записать дважды и взять три последние цифры. 2. Пусть из табл. 2.1 выбраны три группы величин: Коэффициенты СЛАУ Свободные члены c 2,012 –1,347 0,872 1,438 –0,972 3,457 4 + c –2,046 –2,562 1,827 –3,767 6,329 0,899 31 Имея их, надо сформировать СЛАУ 1 2 3 1 2 3 1 2 3 2,012 1,347 0,872 1,438, 3,457 3,028 2,046 2,562, 1,827 3,767 6,329 0,899. x x x x x x x x x           Здесь коэффициент a22 = 3,028 получен сложением 4 и с, равного –0,972. ЛАБОРАТОРНАЯ РАБОТА № 3 РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ИТЕРАЦИОННЫМИ МЕТОДАМИ Цель работы: изучить итерационные методы решения систем линейных алгебраических уравнений – метод простой итерации и метод Зейделя. Состав работы: 1. Изучение алгоритмов итерационных методов решения систем линейных алгебраических уравнений; 2. Пример решения системы линейных уравнений итерационны- ми методами; 3. Выполнение индивидуальных заданий. 3.1. Итерационные методы решения систем линейных алгебраических уравнений Рассмотрим систему линейных алгебраических уравнений n-го порядка 11 1 12 2 1 1 21 1 22 2 2 2 1 1 2 2 ... ... ... . n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b              (3.1) 32 Метод Гаусса предъявляет к решаемой системе единственное требование: она должна быть невырожденной. Для итерационных методов этого требования недостаточно – требуется дополнитель- ное условие, которое называется условием сходимости. Причем имеются два таких условия: 1) необходимое и достаточное; 2) толь- ко достаточное. Однако проверка первого условия часто оказывает- ся более трудоемкой операцией, чем само решение системы. Мы будем пользоваться самым простым, достаточным условием сходимости, одинаковым для обоих методов. Оно называется условием диагонального преобладания и формулируется следую- щим образом: для сходимости итерационного процесса достаточно, чтобы модуль диагонального коэффициента был больше суммы модулей остальных коэффициентов данного уравнения 11 12 1... ,na a a   22 21 23 2... ,na a a a    ……… (3.2) 1 , 1...nn n n na a a    . В отдельных случаях допускается строгое равенство. В математике имеется строгое доказательство того факта, что если диагональные коэффициенты преобладают в каждом уравне- нии системы, то ее определитель отличен от нуля. Обратное утвер- ждение, вообще говоря, неверно. Итак, первым этапом в решении системы уравнений (3.1) итера- ционными методами является проверка условия диагонального преобладания (3.2). Если оно выполняется, переходят ко второму этапу – преобразо- ванию системы к виду, удобному для проведения итераций. Такое преобразование может быть выполнено бесконечно большим чис- лом способов. Наиболее просто это делается путем разрешения пер- вого уравнения относительно x1, второго – относительно x2 и т. д. Получаем: 11 12 1 2 11 11 11 ... ,n n ab ax x x a a a     33 23 22 21 2 1 3 22 22 22 22 ... ,n n a ab ax x x x a a a a     ……… , 11 1 1... n nn n n n nn nn nn ab ax x x a a a      Или 1 1 12 2 1... ,n nx d c x c x    2 2 21 1 23 3 2... ,n nx d c x c x c x     ……… (3.3) 1 1 , 1 1...n n n n n nx d c x c x     , Где / , /i i ii ij ij iid b a c a a  (j ≠ i), cii = 0. (3.4) Система уравнений, приведенная к виду (3.3), называется нор- мальной. На третьем этапе нужно задать начальные значения неизвест- ных (0) (0) (0)1 2, , ... , nx x x , которые называются начальным приближе- нием. Наиболее часто используются два варианта задания начальных значений: 1)  0 0ix  ; 2)  0i ix d , где i = 1, … , n. Однако, если выполняется условие сходимости, итерационный процесс сходится к одним и тем же значениям неизвестных при задании в качестве начального при- ближения любых чисел. На четвертом этапе строится итерационный процесс, состоящий в том, чтобы от любых чисел (0)ix последовательными однотипны- ми шагами (итерациями) приблизиться к истинным значениям xi. Причем на первой итерации (k = 1), используя начальные значения, находят первое приближение (1)ix , на второй итерации (k = 2) – вто- рое приближение (2)ix и т. д. 34 Для вычисления (1)ix в методе простой итерации начальные значения подставляют в правые части всех уравнений нормальной системы (3.3): (1) (0) (0) 1 12 11 2 ... ,n nx d c x c x    (1) (0) (0) (0) 2 21 23 22 1 3 ... ,n nx d c x c x c x     ……… (0) (0)(1) 1 , 11 1...n n n n n nx d c x c x     . Соответственно, для получения второго приближения (2)ix в правые части нормальной системы надо подставить первое приближение: (2) (1) (1) 1 12 11 2 ... ,n nx d c x c x    (2) (1) (1) (1) 2 21 23 22 1 3 ... ,n nx d c x c x c x     ……… (1) (1)(2) 1 , 11 1...n n n n n nx d c x c x     . На k-й итерации (k ≥ 1) находится k-е приближение ( )kix : ( ) ( 1) ( 1) 1 12 11 2 ... , k k k n nx d c x c x      ( ) ( 1) ( 1) ( 1) 2 21 23 22 1 3 ... , k k k k n nx d c x c x c x        ……… (3.5) ( 1) ( 1)( ) 1 , 11 1... k kk n n n n n nx d c x c x       . Метод Зейделя является модификацией метода простой итера- ции и отличается от последнего только способом вычисления при- ближений: ( ) ( 1) ( 1) 1 12 11 2 ... , k k k n nx d c x c x      ( ) ( ) ( 1) ( 1) 2 21 23 22 1 3 ... , k k k k n nx d c x c x c x       ……… (3.6) ( ) ( )( ) 1 , 11 1... k kk n n n n n nx d c x c x     . 35 Это значит, что  1kx находится так же, как и в методе простой итерации, т. е. в правую часть первого уравнения вместо x2, … , xn подставляются значения ( 1) ( 1)2 , ... ,k knx x  , полученные на предыду- щей, (k-1)-ой итерации. При вычислении ( )2kx в правую часть вместо x1 подставляется только что найденное значение ( )1kx , а величины x3, … , xn берутся из предыдущей, (k-1)-ой итерации. Обобщая, можно сказать, что в методе Зейделя значение ( )k ix используются для определения всех последующих величин ( ) ( ) 1 , ... , k k nix x на той же k-ой итерации. При вычислениях по формулам (3.5) и (3.6) возникают следую- щие вопросы: 1) Как ведут себя величины ( )kix по мере увеличения номера итерации k? 2) На какой итерации следует остановить вычислительный про- цесс? 3) Что дает переход от формул (3.5) к формулам (3.6), т. е. от метода простой итерации к методу Зейделя? 4) Что произойдет, если попытаться решить итерационными методами систему, не имеющую диагонального преобладания? 5) Каковы преимущества и недостатки итерационных методов по сравнению с прямыми методами, такими, как метод Гаусса, и како- ва область применения итерационных методов? Ниже даются краткие ответы на эти вопросы. 1. На первых итерациях изменение величин ( )kix носит, как правило, хаотичный характер. Однако постепенно их значения ста- билизируются: сначала устанавливается правильный знак числа, затем появляется и начинает повторяться первая правильная цифра, затем вторая и т. д. Это свидетельствует о том, что итерационный процесс сходится. Причем сходится независимо от того, с какого начального приближения он стартовал. 2. Процесс вычислений по формулам (3.5) и (3.6) является в принципе бесконечным. И остановить его надо в тот момент, когда будет получено решение с наперед заданной точностью 36 (m верных значащих цифр). Исследования показывают, что необхо- димое для этого число итераций зависит не только от требуемой точности результата, но и от свойств самой системы (большей или меньшей степени выраженности ее диагонального преоблада- ния, порядка), а также величины свободных членов и начальных значений неизвестных. Однако формулы и методики определения минимального числа итераций сложны, особенно в случае метода Зейделя. Поэтому на практике лучше пользоваться более простыми приемами. Если система невысокого порядка решается с использованием калькулятора с последовательной записью на бумаге всех значений ( )kxi , то имеется возможность в процессе решения наблюдать за их поведением. При этом легко обнаруживается тот момент, когда все неизвестные в двух соседних итерациях совпадают друг с другом до m +1 верных значащих цифр. Тогда, округлив ( )kxi до m цифр, гарантировано получают результат требуемой точности. При счете на ПК следить визуально за ходом вычислений прак- тически невозможно. Поэтому на каждой итерации надо выполнять проверку ( ) ( 1) ( ) k k i i k i x x x    , (3.7) где  – некоторое малое, наперед заданное число. Вычислительный процесс прекращается, если условие (3.7) выполняется для всех i от 1 до n. Это условие не срабатывает в одном случае, а именно, когда ( )k ix окажется равным нулю. Чтобы предупредить этот случай, представим (3.7) в виде ( ) ( 1) ( )k k k i i ix x x    , а затем добавим к правой части малую величину  ( ) ( 1) ( )k k k i i ix x x      . (3.8) 37 Это условие будет «работать» при любом ( )kix . В частности, при ( ) 0kix  оно вырождается в следующее: ( 1)k ix    . Что касается величины , то анализ показывает: для получения результата с m верными значащими цифрами надо принять  не более, чем 5·10-(m+2). 3. Интуиция подсказывает, что немедленный ввод в вычисления найденного на этой же итерации значения неизвестного, как это предусматривают формулы (3.6), должен привести к ускорению сходимости. И практика подтверждает, что метод Зейделя требует, как правило, меньшего числа итераций для получения результата той же точности. Однако строгого доказательства этого факта пока нет. 4. Если попытаться решать итерационными методами систему, не имеющую диагонального преобладания, то вычислительный процесс будет, как правило, расходящимся. Это значит, что с уве- личением номера итерации значения неизвестных не стабилизиру- ются, а неограниченно возрастают со знаком плюс, или минус, или с переменой знака. Наконец, наступает момент переполнения памяти машины и происходит остановка выполнения программы. Однако такой исход вовсе не обязателен. Дело в том, что условие диагонального преобладания является только достаточным. Кроме того, области сходимости методов простой итерации и Зейделя пе- ресекаются, но полностью не совпадают. Поэтому при рассмотре- нии конкретной системы могут встретиться четыре случая: 1) система не решается обоими итерационными методами; 2) система решается только методом простой итерации; 3) система решается только методом Зейделя; 4) система решается обоими методами. Как мы уже знаем, условие диагонального преобладания гаран- тирует четвертый случай. 5. Итак, область применения итерационных методов – решение, главным образом, систем с диагональным преобладанием. Это недостаток итерационных методов, причем довольно суще- ственный. Хотя в расчетной практике такие системы встречаются довольно часто. Однако эти методы имеют и преимущества. 38 Например, они в отличие от прямых методов дают результаты за- данной точности. Поэтому итерационные методы можно использо- вать для уточнения решения, полученного прямыми методами (при этом значения неизвестных, найденные прямым методом, берутся в качестве начального приближения для итерационного уточнения). В практических задачах возможно появление такой ситуации, когда система имеет настолько высокий порядок, что ее невозмож- но расположить в оперативной памяти ПК, но уравнения очень про- стые. В этом случае их формирование можно непосредственно включить в итерационный процесс решения, что при применении прямых методов сделать невозможно. 3.2. Пример решения системы линейных алгебраических уравнений итерационными методами Снова рассмотрим систему линейных уравнений третьего порядка 1 2 3 1 2 3 1 2 3 25,2 4,7 1,4 46,8, 4,1 39,4 2,9 25,6, 1,3 3,1 16,8 34,3 , x x x x x x x x x            (3.9) которая уже была решена в лабораторной работе № 2 методом Гаус- са. В настоящей работе найдем решение этой системы итерационны- ми методами – методом простой итерации и методом Зейделя. Решение. 1. Проверяем условие вырожденности системы. Так как опреде- литель матрицы не равен 0, то решение системы существует и един- ственно. 2. Проверим условие диагонального преобладания. 25,2 4,7 1,4   ; 39,4 4,1 2,9  ; 1,33,18,16  . Так как условие выполнятся, то для поиска решения можно ис- пользовать оба метода: простой итерации и Зейделя. 3. Преобразование заданной системы к нормальному виду. Решая первое уравнение относительно x1, второе – относительно x2, а третье – относительно x3, получаем 39 x1 = 1,85714 + 0,186508x2 – 0,0555556x3, x2 = –0,649746 – 0,104061x1 – 0,0736041x3, (3.10) x3 = 2,04167 + 0,0773810x1 – 0,184524x2. 4. Задание начального приближения. Мы уже имеем решение этой системы методом Гаусса. И его можно бы использовать в качестве начального приближения для итерационного уточнения решения. Однако, чтобы показать возможности итерационных методов, примем в качестве начальных значений числа, которые отличаются от искомых величин не только абсолютным значением, но и знаком: (0)1 3,x   (0)2 2,x  (0)3 1x  . 5. Проведение итераций. Вначале применим метод простой итерации. Подставив начальные значения в правую часть нормальной системы (3.10), получаем первое приближение (1) 1 1,85714 0,186508 2 0,0555556 1 2,17460,x       (1) 2 0,649746 0,104061 ( 3) 0,0736041 1 0,411168,x          (1) 3 2,04167 0,0773810 ( 3) 0,184524 2 1,44048x        . Полученные значения вновь подставляем в правую часть нормальной системы и находим второе приближение (2) 1 1,85714 0,186508 ( 0,411168) 0,0555556 1,44048 1,70043,x        (2) 2 0,649746 0,104061 2,17460 0,0736041 1,44048 0,982062,x         (2) 3 2,04167 0,0773810 2,17460 0,184524 ( 0,411168) 2,28581.x        Эти и последующие приближения неизвестных сведем в табл. 3.1, где величина k обозначает номер итерации или номер приближения. 40 Таблица 3.1 Решение системы (3.9) методом простой итерации k ( )1 kx ( )2 kx ( )3 kx 0 –3 2 1 1 2,17460 –0,411168 1,44048 2 1,70043 –0,982062 2,28581 3 1,54699 –0,994939 2,35446 4 1,54078 –0,984025 2,34396 5 1,54334 –0,982680 2,34247 6 1,54373 –0,982763 2,34242 7 1,54372 –0,982800 2,34246 8 1,54371 –0,982802 2,34247 Анализ этой таблицы показывает, как формируется решение в процесс сходимости: правильные знаки неизвестных устанавли- ваются уже на первой итерации, первая правильная цифра – на вто- рой итерации, вторая – на четвертой итерации и т. д. Процесс вычислений остановлен в тот момент, когда у всех неизвестных повторились первые пять значащих цифр. Округлив эти значения до четырех цифр, получаем x1 = 1,5437, x2 = –0,9828, x3 = 2,3425. Найденные значения, вообще говоря, подлежат проверке подста- новкой их во все уравнения заданной системы, как это делалось в лабораторной работе № 2. Однако факт совпадения полученных значений с решением, найденным методом Гаусса, тоже свидетель- ствуют о правильности результатов. Решим теперь исходную систему по методу Зейделя, вычисляя приближения по формулам (3.6). Имеем первое приближение: (1) (0) (0) 1 12 131 2 3 1,85714 0,186508 2 0,0555556 1 2,17460,x d c x c x         (1) (1) (0) 2 21 232 1 3 0,649746 0,104061 2,17460 0,0736041 1 0,949641, x d c x c x            (1) (1) (1) 3 31 323 1 2 2,04167 0,0773810 2,17460 0,184524 ( 0,949641) 2,38517. x d c x c x           41 Второе приближение (2) (1) (1) 1 12 131 2 3 1,85714 0,186508 ( 0,949641) 0,0555556 2,38517 1,54572, x d c x c x           (2) (2) (1) 2 21 232 1 3 0,649746 0,104061 1,54752 0,0736041 2,38517 0,986341, x d c x c x            (2) (2) (2) 3 31 323 1 2 2,04167 0,0773810 1,54752 0,184524 ( 0,986341) 2,34342. x d c x c x           За поведением решения можно проследить по табл. 3.2. Таблица 3.2 Решение системы (3.9) методом Зейделя k ( )1 kx ( )2 kx ( )3 kx 0 –3 2 1 1 2,17460 –0,949641 2,38517 2 1,54752 –0,986341 2,34342 3 1,54299 –0,982796 2,34247 4 1,54371 –0,982797 2,34247 5 1,54371 –0,982801 2,34247 Как и следовало ожидать, итерационный процесс сошелся быст- рее, чем в методе простой итерации. 3.3. Выполнение индивидуальных заданий Составить систему линейных алгебраических уравнений, выбрав коэффициенты и свободные члены из табл. 3.3 согласно шифру. Решить эту систему методом простой итерации и методом Зейделя с точностью E = 10–4. Сделать проверку полученных результатов. Решить задание, используя систему Mathematica, и сравнить полу- ченные результаты. 42 Таблица 3.3. Исходные данные к заданию Пер- вая цифра шифра Коэффициенты системы уравнений Вто-рая цифра шиф- ра Свобод- ные члены bi Третья цифра шифра Вели- чина с ai1 ai2 ai3 1 2 3 4 5 6 7 8 0 14,34 – c –4,157 –3,306 1,588 12,750 4,530 –2,164 2,949 9,760 0 0,670 –4,174 0,428 0 0,513 1 11,230 4,679 –8,417 1,478 15,208 + c 0,784 3,245 1,039 10,732 1 1,767 3,134 1,987 1 0,317 2 9,214 5,298 4,297 3,047 7,319 1,040 2,046 1,395 7,943 2 4,278 1,409 0,903 + c 2 1,783 3 23,729 8,860 3,495 10,395 – c 14,395 2,845 1,010 4,045 13,230 3 2,170 4,194 3,294 3 –3,134 4 5,204 1,054 +c 0,371 3,956 20,423 5,320 0,030 2,569 17,456 4 4,276 –3,167 1,328 4 0,451 5 –13,057 1,390 5,281 2,239 19,432 4,130 + c 1,945 4,853 19,634 5 13,045 2,056 –1,940 5 0,574 6 11,945 2,478 3,287 1,035 15,374 7,813 3,275 4,230 +c 13,743 6 2,874 4,934 5,397 6 0,323 7 14,023 2,056 1,392 0,456 12,385 4,439 1,493 +c 5,179 13,845 7 6,905 2,504 4,05 7 0,124 8 9,394 9,010 2,156 1,034 15,345 7,032 – c 0,638 4,945 19,045 8 3,905 5,234 1,043 + c 8 –0,394 9 13,056 –3,456 –5,170 0,045 – c 7,845 –6,134 –0,198 –0,345 15,390 9 3,560 – c 2,157 10,150 9 –1,489 43 ЛАБОРАТОРНАЯ РАБОТА № 4. НАХОЖДЕНИЕ НАИБОЛЬШЕГО ПО МОДУЛЮ СОБСТВЕННОГО ЗНАЧЕНИЯ И СООТВЕТСТВУЮЩЕГО ЕМУ СОБСТВЕННОГО ВЕКТОРА КВАДРАТНОЙ МАТРИЦЫ Цель работы: изучить степенной итерационный метод и его модификации для определения наибольшего по модулю собствен- ного значения 1 и соответствующего ему собственного вектора 1x квадратной матрицы. Состав работы: 1. Изучение алгоритма степенного итерационного метода. 2. Пример определения 1 и 1x степенным методом; 3. Недостатки степенного итерационного метода; 4. Модификации степенного итерационного метода; 5. Пример определения 1 и 1x модифицированным степенным методом; 6. Поиск собственных значений и собственных векторов матри- цы в системе Mathematica; 7. Выполнение индивидуальных заданий. 4.1. Нахождение 1 и 1x степенным методом Пусть дана квадратная матрица A произвольного порядка n. Требуется решить для нее с заданной точностью частичную про- блему собственных значений по определению 1 и 1x , т. е. такого числа 1 и такого вектора 1x , которые удовлетворяют соотношению 1 1 1x x A   (4.1) (как известно, аналогичным соотношениям удовлетворяют второе и последующие собственные значения: 2 2 2x x A   , …). Для решения этой задачи степенным итерационным методом нужно задаться начальным вектором (0)y порядка n , имеющим 44 произвольные, но ненулевые компоненты. Затем строится итераци- онный процесс ( ) ( 1)k ky y  A  (k = 1, 2, … – номера итераций), в результате которого образуется последовательность векторов (0)y  , (1) (0)y y A  , (2) (1)y y A  , … . В ходе вычисления этих векторов можно находить все более уточняющееся значение 1 по формуле ( ) ( ) 1 ( 1) k k s k s y y    (k = 1, 2, …), (4.2) где ( )ksy – s-я (любая, s = 1, … , n) компонента вектора ( )ky на k-й итерации; ( 1)k sy  – соответствующая компонента вектора ( 1)ky  на преды- дущей итерации. k-е приближение для первого собственного вектора можно брать как ( ) ( )k k ix y  (k = 1, 2, …). (4.3) Для доказательства формул (4.2) и (4.3) разложим начальный вектор (0)y по собственным векторам матрицы A (0) 1 1 2 2 n ny a x a x a x       , (4.4) где ai (i = 1, … , n) – некоторые числа, причем a1 ≠ 0. Тогда (1) (0) 1 1 2 2 n ny y a x a x a x    A A A A     . 45 Исходя из определения собственного значения и собственного вектора 1 1 1x x A   , 2 2 2x x A   , … , n n nx x A   . (4.5) Поэтому (1) 1 1 1 2 2 2 n n ny a x a x a x          . Переходя ко второму вектору, получаем (2) (1) 1 1 1 2 2 2 n n ny y a x a x a x       A A A A     или (2) 2 2 2 1 1 1 2 2 2 n n ny a x a x a x          . (4.6) Теперь очевидно, что вектор ( )ky (k = 1, 2, …) содержит в своем разложении собственные значения, возведенные в k-ю степень (отсюда и название рассматриваемого метода – степенной). Но по- скольку 1 – наибольшее из всех n собственных значений матрицы, то доля первого слагаемого в разложениях (4.5), (4.6) и т. д. с уве- личением номера итерации будет все более возрастать. Это дает основание в записанных выше и последующих разложениях отбро- сить второе, третье, … , n-е слагаемое, записав приближенные равенства (0) 1 1y a x  , (1) 1 1 1y a x   , (4.7) (2) 2 1 1 1y a x   , ……… Такие же зависимости справедливы и для отдельных компонен- тов входящих сюда векторов (0) 1 1s sy a x , (1) 1 1 1s sy a x  , (2) 2 1 1 1s sy a x  , ……… , где s может принимать любое значение среди чисел 1, … , n. 46 Разделив второе равенство на первое, получим 1 в первом при- ближении (1) (1) 1 (0) s s y y   . Чтобы получить второе приближение для 1, нужно разделить третье равенство на второе: (2) (2) 1 (1) s s y y   . На любой k- й итерации имеем формулу (4.2). Справедливость формулы (4.3) вытекает из приближенных равенств (4.7). При этом следует вспомнить, что собственный вектор остается таким же, будучи умноженным на любое число, отличное от нуля. Поскольку a1 ≠ 0,то на первой итерации за 1x можно взять вектор 1 1 1a x  или 1y , на второй – 21 1 1a x  или (2)y и т. д. Конечно, проводя вычисления по формулам (4.2) и (4.3), мы по- лучим для 1 и 1x на первых итерациях самые немыслимые значе- ния. Однако постепенно искомые величины будут стабилизировать- ся, пока не сойдутся с заданной точностью к своим истинным зна- чениям, что можно проверить подстановкой в соотношение (4.1). Момент окончания счета определяется так же, как и в итерационных методах решения систем линейных алгебраических уравнений. 4.2. Пример определения 1 и 1x степенным методом Пусть дана матрица третьего порядка 25,2 4,7 1,4 4,1 39,4 2,9 1,3 3,1 16,8        A . (4.8) Требуется определить 1 и 1x с точностью  = 10-3, используя степенной итерационный метод. 47 Решение. Зададимся начальным вектором третьего порядка  0 5 2 1 Тy   и вычислим (1) (0) 25,2 4,7 1,4 5 136,8 4,1 39,4 2,9 2 55,4 1,3 3,1 16,8 1 4,1 y y                                A  . Примем s = 2. Тогда первое приближение для 1 определится как (1) (1) (0) 1 2 2/ 55,4 / ( 2) 27,7y y      . Соответствующий собственный вектор  (1) (1)1 136,8 55,4 4,1 Тx y   . Переходим ко второй итерации: (2) (1) 25,2 4,7 1,4 136,8 3713,48 4,1 39,4 2,9 55,4 1609,99 1,3 3,1 16,8 4,1 280,7 y y                                 A  , (2) (2) (1) 1 2 2/ 1609,99 / ( 55,4) 29,0612y y      ,  (2) (2)1 3 713,48 1609,99 280,7 Тx y    . Продолжая итерационный процесс, сведем эти и последующие результаты в табл. 4.1. Таблица 4.1 Нахождение 1 и 1x для матрицы (4.8) степенным итерационным методом k ( ) 1 ky ( )2 ky ( )3 ky  1k  = ( ) ( 1)2 2/ k ky y   k  =    1 1 1 k k   0 5 –2 1 – – 1 1,368·102 –5,54·101 4,1·100 27,7 – 2 3,71348·103 –1,60999·103 –2,807·102 29,0612 1,3612 48 Окончание табл. 4.1 1 2 3 4 5 6 3 1,00754·105 –4,90224·104 –1,45343·104 30,4489 1,3877 … … … … … … 26 2,29435·104 0 –6,81085·1040 –1,11122·1040 38,4912 0,0012 27 8,82730·104 1 –2,62163·1042 –4,27648·1041 38,4920 0,0008 Вычисления прекращены на 27-й итерации, так как (27) меньше  = 10-3: (27) 1 38,4920  , (27) 41 42 41 1 8,82730 10 2,62163 10 4,27648 10 Т x          . Проверим эти значения подстановкой в равенство (4.1): 41 43 (27) 42 44 1 41 43 8,82730 10 3,39678 1025,2 4,7 1,4 4,1 39,4 2,9 2,62163 10 1,00913 10 1,3 3,1 16,8 4,27648 10 1,64591 10 x                                       A  , 41 43 (27) (27) 42 44 1 1 41 43 8,82730 10 3,39780 10 38,4920 2,62163 10 1,00912 10 4,27648 10 1,64610 10 x                                . Как видно, соответствующие компоненты левой и правой частей равенства (4.1) практически совпадают друг с другом. Таким обра- зом, полученный результат следующий: 1 = 38,492; 41 42 411 8,8273 10 2,62163 10 4,27648 10 Т x          . Однако собственным вектором, представленным таким образом, пользоваться неудобно. Чаще всего его тем или иным способом нормируют. Наиболее употребительны два способа нормирования: 49 1) делением 1x на его наибольшую по модулю компоненту (в данном случае таковой является вторая компонента – 2,62163·1042). Тогда в новом (нормированном) векторе соответ- ствующая компонента станет равной +1, а остальные по мо дулю не будут превосходить единицу. Нормированный таким спо- собом собственный вектор обозначим одной звездочкой:  *1 0,3367 1 0,1631 Tx   ; 2) делением 1x на его длину (модуль). Полученный таким спосо- бом вектор называется нормированным к единичной длине. Мы будем обозначать его двумя звездочками. В нашей задаче (27) 41 2 42 2 41 2 42(8,82730 10 ) ( 2,62163 10 ) ( 4,27648 10 ) 2,79912 10x           и  **1 0,3154 0,9366 0,1528 Tx    . Длина этого вектора равна 1. Векторы 1x , *1x , **1x отличаются друг от друга на некоторый множитель, т. е. имеют разную длину, но одинаковое направление, поэтому каждый из них является собственным вектором заданной матрицы (4.8). 4.3. Недостатки степенного итерационного метода Степенной итерационный метод обладает рядом существенных недостатков: 1) если 1 1  (как это имеет место в нашем примере), то с уве- личением номера итерации k компоненты вектора ( )ky могут воз- расти настолько, что превысят верхнюю границу области допусти- мых значений для вещественных чисел и произойдет аварийная остановка программы при ее выполнении на ЭВМ; 2) если 1 1  , то с увеличением номера итерации k компоненты вектора ( )ky могут уменьшиться настолько, что они окажутся мень- 50 ше нижней границы области допустимых значений для веществен- ных чисел и будут заменены машинным нулем, что означает потерю информации; 3) В ходе итерационного процесса компоненты вектора ( )ky могут менять знак (как это случилось в нашем примере с третьей компонентой при переходе от первого приближения ко второму, см. табл. 4.1). В процессе перемены знака компонента может пройти через нуль, и тогда при делении, как это предусматривает формула (4.2), произойдет аварийная остановка программы. Эта же причина вынуждает отказаться от задания нулевых компонент в начальном векторе (0)y ; 4) Сходимость процесса можно отслеживать не только по 1. Между тем, компоненты собственного вектора могут иметь более медленную сходимость, что может быть установлено только при проверке результатов. Проще всего устраняется третий недостаток. Для этого ( )1k надо вычислять как среднее из тех m частных ( ) ( 1)/k ks sy y  , у которых ( 1) 0ksy   , т. е. ( ) ( ) 1 ( 1) 1 1 knk s k s s y m y     . Здесь m ≤ n, т. е. слагаемые, которым отвечает ( 1) 0ksy   , про- пускаются. Теперь к заданию компонент вектора (0)y может использоваться менее жесткий подход, однако все они не могут быть одновременно равными нулю. 4.4. Модификации степенного итерационного метода Существенные недостатки, присущие степенному итерационно- му методу, явились предпосылкой исследования новых путей его совершенствования. Так появились модифицированные степенные методы. 51 Самая простая модификация, ликвидирующая, однако, сразу все четыре недостатки степенного итерационного метода (п. 4.3), состоит в том, что последовательность векторов строится по формулам (0)y , (0) (0)(0) 1 s x y y   , (1) (0)y x A  , (1) (1)(1) 1 s x y y   , (2) (1)y x A  , (2) (2)(2) 1 s x y y   , ……… . Здесь ( )ksy – наибольшая по модулю компонента вектора ( )ky . Это значит, что на каждой итерации выполняется нормирование вектора первым способом, т. е. делением его на наибольшую по мо- дулю компоненту. И этот нормированный вектор используется для нахождения следующего вектора ( 1)ky  . При таком построении итерационного процесса формула (4.2) примет вид ( ) ( ) ( 1) 1 / k k k s sy x   , (k = 1, 2, …), (4.11) где ( )ksy – наибольшая по модулю компонента вектора ( )ky , а ( 1)ksx  – соответствующая компонента вектора ( 1)kx  . На одной-двух первых итерациях номер s компоненты, имеющей наибольший модуль, может меняться. Дальше s устанавливается и сохраняет свое значение, после чего на всех итерациях ( 1)ksx  будет равен единице и формула (4.11) упрощается ( ) ( ) 1 k k sy  , (k = 1, 2, …). (4.12) Мы будем пользоваться формулой (4.12) и на первых итерациях, т.к. полученные на них значения 1 будут далекими от истинной величины независимо от того, по какой формуле их вычислять. 52 В качестве k-го приближения первого собственного вектора мож- но взять ( )kx . Это будет нормированный первым способом вектор ( ) ( ) 1 k kx x  , (k = 1, 2, …). (4.13) Можно построить и другой модифицированный степенной метод, нормируя векторы ( )ky вторым способом (делением на длинну вектора ( )ky . В этом случае строится последователь- ность векторов: (0)y , (0) (0) (0) 1x y y   , (1) (0)y A x  , (1) (1) (1) 1x y y   , (4.14) (2) (1)y A x  , (2) (2) (2) 1x y y   , ……… . При этом k-е приближение для первого собственного значения определяется по формуле ( ) ( ) ( 1) 1 ( , ) k k ky x     , (k = 1, 2, …), что означает скалярное произведение двух векторов, равное, как известно, сумме произведений их соответствующих компонент, т. е. ( ) ( ) ( 1) 1 1 nk k k i i i y x      , (k = 1, 2, …). (4.15) k-е приближение для первого собственного вектора, нормированно- го к единице, находится как **( ) ( ) 1 k kx x  , (k = 1, 2, …). (4.16) 53 4.5. Пример определения 1 и 1x модифицированным степенным методом Рассмотрим снова матрицу (4.8) и вычислим для нее 1 и 1x модифицированным степенным методом с точностью  = 10-3. Зададим тот же начальный вектор  0 5 2 1 Тy   . Его наиболь- шая по модулю компонента (0)1 5y  и имеет номер s = 1. Тогда начальный нормированный вектор будет  (0) 1 0,4 0,2 Тx   . Выполним первую итерацию (1) (0) 25,2 4,7 1,4 1 27,36 4,1 39,4 2,9 0,4 11,08 1,3 3,1 16,8 0,2 0,82 y x                                A  , (1) (1) 1 1y  = 27,36 (наибольшая компонента, s = 1). Нормированный вектор   *(1)(1) 11 0,405 0,03 Тx x    . Перейдем ко второй итерации (2) (1) 25,2 4,7 1,4 1 27,1453 4,1 39,4 2,9 0,405 11,7689 1,3 3,1 16,8 0,003 2,0519 y x                                 A  , (2) (2) 1 1y  = 27,1453. Нормированный вектор   *(2)(2) 11 0,4336 0,0756 Тx x     . Эти и все последующие значения на итерациях сведены в табл. 4.2. 54 Таблица 4.2 Нахождение 1 и 1x для матрицы (4.8) модифицированным степенным методом k ( ) 1 ky ( )2 ky ( )3 ky  1k    max ky  k     1k k   ( ) 1 kx ( )2 kx ( )3 kx 0 5 –2 1 1 –0,4 0,2 1 27,36 –11,08 0,82 27,36 1 –0,405 0,03 2 27,1453 –11,7689 –2,0519 27,1453 0,2147 1 –0,4336 –0,0756 3 27,1319 –13,2012 –3,9139 27,1319 0,0134 1 –0,4866 –0,1443 … … … … … … … … … 6 28,0866 –23,4296 –7,2375 28,0866 0,4871 1 –0,8342 –0,2577 7 28,7599 –29,5144 –8,2151 28,7599 0,6733 –0.9744 1 0,2783 8 –28,8662 36,2120 9,0429 36,2120 7,4521 –0,7971 1 0,2497 … … … … … … … … … 26 –12,9664 38,4912 6,2800 38,4912 0,0012 –0,3369 1 0,1632 27 –12,9606 38,4920 6,2789 38,4920 0,0008 –0,3367 1 0,1631 Следует отметить, что наибольшая компонента на 7-й итерации поменяла номер (было s = 1, стало s = 2). Остановив вычисления при (k) ≤  (это 27-я итерация), получаем ответ, аналогичный как и при решении предыдущим способом (27) 1 2y  = 38,492;  * (27)1 0,3367 1 0,1631 Тx x    . 4.6. Поиск собственных значений и собственных векторов матрицы в системе Mathematica Собственные значения и собственные вектора матриц можно найти, используя функции Eigenvalues и Eigenvectors. Для того, чтобы определить характеристический полином матрицы, ис- пользуется функция CharacteristicPolynomial. Реализация примера из п. 4.2 в системе Mathematica показана на рисунке. 55 Функции Eigenvalues и Eigenvectors возвращают все собственные значения и соответствующие им соб- ственные вектора. Значение рассчитанного в системе Mathematica наибольшего по модулю значения  равно 38,4937, что практически совпадает с найденным с по- мощью итерационных мето- дов. Значения компонент собственного вектора также практически совпадают (см. п. 4.2 при использование нормировки на длину век- тора). 4.7. Выполнение индивидуальных заданий Составить матрицу, выбрав коэффициенты из табл. 4.3 согласно шифру. Используя степенной итерационный метод и его модифика- цию, определить наибольшее по модулю собственное значение  (точность  = 10-3) и соответствующий ему собственный вектор. Провести нормировку собственного вектора двумя способами. Сделать проверку полученных результатов. Решить задание, используя систему Mathematica. Таблица 4.3 Исходные данные к заданию Пер- вая цифра шифра Коэффициенты системы уравнений Вто- рая цифра шифра Сво- бод- ные члены bi Третья цифра шифра Вели- чина с ai1 ai2 ai3 1 2 3 4 5 6 7 8 0 1,034 15,345 + с 2,056 2,874 – d 0,456 5,204 2,504 11,945 13,056 0 0,124 0 1,783 1 4,157 12,750 – d 8,860 0,371 4,218 1,390 0,030 3,275 + с 1,390 1 0,574 1 0,030 56 Окончание табл. 4.3 1 2 3 4 5 6 7 8 2 1,588 + с 14,395 1,390 4,945 4,276 3,495 13,045 1,035 5,179 + d 2 1,040 2 –1,940 3 14,023 2,478 2,239 + с 5,204 5,281 + d 5,32 2,874 9,010 17,456 3 0,323 3 –0,394 4 13,057 5,204 + d –3,456 15,345 5,397 3,287 7,813 1,430 + с 4,130 4 –0,198 4 1,588 5 0,371 + с 12,385 2,949 4,853 13,729 – d 4,045 2,845 –0,345 1,034 5 0,638 5 0,317 6 4,230 9,394 2,156 – d 3,956 + с 15,374 –6,134 1,493 3,245 –4,157 6 0,451 6 1,054 7 8,860 11,945 19,045 19,432 7,032 + с 2,056 0,456 7,319 + d 19,634 7 –1,489 7 1,395 8 2,056 2,569 + с 12,750 4,439 7,845 13,230 0,045 – d 13,845 3,275 8 1,010 8 –3,134 9 2,046 5,170 4,194 + с 13,056 4,934 20,423 1,392 1,945 – d 4,679 9 0,513 9 0,039 57 ЛАБОРАТОРНАЯ РАБОТА № 5 РЕШЕНИЕ ЗАДАЧИ КОШИ МЕТОДАМИ РУНГЕ-КУТТА Цель работы: изучить методы Рунге-Кутта решения одного дифференциального уравнения первого порядка, систем дифферен- циальных уравнений первого порядка, а также уравнений второго и старших порядков. Состав работы: 1. Общие сведения о дифференциальных уравнениях; 2. Решение задачи Коши для одного дифференциального уравне- ния первого порядка; 3. Пример решения задачи Коши для одного дифференциального уравнения первого порядка; 4. Применение методов Рунге-Кутта к решению систем диффе- ренциальных уравнений первого порядка; 5. Пример решения задачи Коши для системы двух линейных дифференциальных уравнений первого порядка; 6. Применение методов Рунге-Кутта к решению дифференциаль- ных уравнений второго и старшего порядков; 7. Решение задачи Коши в системе Mathematica. 8. Выполнение индивидуальных заданий. 5.1. Общие сведения о дифференциальных уравнениях Обыкновенными дифференциальными уравнениями называются уравнения, содержащие одну независимую переменную, искомую функцию y = y(x) и одну или несколько ее производных. Например, F(x, y, y', … ,y(n)) = 0. Здесь x – независимая переменная, y(i), i = 1, 2, … n – i-я производная неизвестной функции y(x); F – некоторая функция от переменных x, y, y', y'', … . Порядок дифференциального уравнения определяется наивыс- шим порядком n входящей в уравнение производной. Если в урав- нении старшую производную выразить через производные более низкого порядка и независимую переменную, то получаются урав- нения, разрешенные относительно старшей производной. Напри- мер, y''' = f(x, y, y',y''). Дифференциальные уравнения бывают линейные и нелинейные. Линейным называется уравнение, в котором неизвестная функция 58 и ее производные входят в первой степени. Например, y'' = 2xy + x2. Остальные уравнения называются нелинейными. Решением дифференциального уравнения называется произ- вольная функция (x), которая при подстановке в уравнение пре- вращает его в тождество. Общее решение обыкновенного дифференциального уравнения n-го порядка, как известно, содержит n произвольных постоянных c1, c2, … cn. Частное решение получается из общего путем задания определенных значений произвольным постоянным ci, i = 1, 2, … n. Для выделения частного решения из общего нужно задать столько дополнительных условий, сколько произвольных постоянных со- держится в общем решении, т. е. каков порядок дифференциального уравнения. Для системы n обыкновенных дифференциальных уравнений первого порядка в качестве дополнительных условий задаются зна- чения искомых функций и их производные при некоторых x. В зависимости от способа задания дополнительных условий для получения частного решения дифференциального уравнения разли- чают два типа задач: задача Коши и краевые задачи. Если дополнительные условия задаются в одной точке x, то зада- ча называется задачей Коши, а дополнительные условия называют- ся начальными условиями. Точка х = x0, в которой они задаются, – начальной точкой. Для уравнения первого порядка задается одно начальное условие y(x0) = y0, т. е. начальное значение искомой функции. Для уравнения второго порядка задаются два начальных условия y(x0) = y0 и y'(x0) = y0', т. е. начальное значение искомой функции и ее первой производной. В общем случае для уравнения n-го порядка задаются: y(x0) = y0, y'(x0) = y0', … y(n-1)(x0) = y0(n-1). Обычно задача Коши связана с изучением процессов, протекаю- щих во времени (колебания строительных конструкций, процесс распределения тепла и т. д.). Поэтому аргумент в таких задачах (время) обозначается через t, а производные записываются с использованием точки. Если дополнительные условия задаются в нескольких точках, то есть при разных х, то такая задача называется краевой, а допол- нительные условия – краевыми или граничными условиями. Обычно граничные условия задаются в точках х = а и х = b, являю- щихся границами области решения дифференциального уравнения. 59 5.2. Решение задачи Коши для одного дифференциального уравнения первого порядка Пусть имеется задача Коши  ,y f t y , y(t0) = y0. Требуется решить ее численно, т. е. найти значения y1, y2, … при значениях аргумента t1, t2, … , отстоящих друг от друга на шаг ∆t. Это означает, что t1 = t0 + ∆t; t2 = t1 + ∆t, … . Все многообразие методов решения этой задачи делится на две группы – одношаговые и многошаговые. В одношаговых методах величина yi определяется по информации, имеющейся в предыду- щем узле i – 1. Соответственно, в многошаговых методах использу- ется информация в k предыдущих узлах (k ≥ 2). Точность многоша- говых методов, как правило, выше, но реализовать их сложнее. Поэтому в расчетной практике используют главным образом одно- шаговые методы. Самым простым и исторически первым является метод лома- ных Леонардо Эйлера. Вместе с тем, это и самый неточный метод. В нем очередное значение искомой функции определяется по формуле yi = yi–1 + ki–1∆t,(5.1) где ki–1 – тангенс угла наклона касательной к графику искомой функции в точке (ti–1, yi–1). Чтобы его определить, надо величины ti–1 и yi–1 подставить в правую часть решаемого уравнения соответ- ственно вместо t и y, т. е. ki–1 = f(ti–1, yi–1). Геометрическая интерпретация метода показана на рис. 5.1. К интегральной кривой, проходящей через точку A с координа- тами ti–1 и yi–1, проведена касательная AB. Тангенс угла наклона ее к горизонту i–1 обозначен через ki–1. 60 Очередное значение искомой функции yi складывается из преды- дущего значения yi–1 и приращения BC, которое можно определить следую-щим образом |BC| = |AC|tg(i–1) = ki–1∆t. Из рисунка видно, что найденное значение yi отличается от истинного значения на величину BD, являющуюся ошибкой ме- тода, и наряду с ошибкой округления делающую результат прибли- женным. Ошибки способны накапливаться, снижая точность реше- ния с каждым шагом. Единственный способ повысить точность решения в методе Эйлера – уменьшить шаг. Именно это обстоятельство побудило немецких математиков К. Рунге и В. Кутта разработать целую гам- му одношаговых методов, отличающихся порядком. Сюда входят: метод первого порядка (полностью совпадает с методом Эйлера), метод второго порядка и т. д. Чем выше порядок метода, тем точнее результат (при одном и том же шаге), но тем больше арифметиче- ских действий надо выполнить. В методе Рунга-Кутта второго порядка вначале делается пробный шаг в узел i, как и методе первого порядка: ki–1 = f(ti–1, yi–1), yi* = yi–1 + ki–1∆t. t A D B С ttiti-1 yi-1 yi-1 i ki-1t Интегральная кривая yi Ошибка метода ki-1=tg(i Рис. 5.1. Геометрическая интерпретация метода Эйлера 61 Затем находят тангенс угла i* наклона касательной к интеграль- ной кривой в точке (ti, yi*) ki* = f(ti, yi*) и определяют более точное значение искомой функции yi = yi–1 + 1/2·(ki–1 + ki*)·∆t(5.2) (касательная из точки A проводится под углом, тангенс которого равен полусумме тангенсов ki-1 и ki*). В настоящее время широкое применение в расчетной практике находит также метод Рунге-Кутта четвертого порядка. 5.3. Пример решения задачи Коши для одного дифференциального уравнения первого порядка Решить методами первого и второго порядка задачу Коши 2y y t , y(1) = 0,5(5.3) на отрезке [1, 2] в узлах с шагом ∆t = 0,2. Эта задача имеет точное аналитическое решение y(t) = 0,5t2, которое позже используем для сравнения с приближенными резуль- татами. Решим эту задачу методом Эйлера. Имеем: 1. Для первого узла (i = 1) 0 0 0 0,52 2 1 1 yk t     ; 1 0 0 0,5 1 0, 2 0,7y y k t       ; 2. Для второго узла (i = 2) 1 1 1 0,72 2 1,16667 1,2 yk t     ; 2 1 1 0, 7 1,16667 0, 2 0,933333y y k t       . Аналогично проводятся вычисления для остальных узлов. Полученные результаты показаны в табл. 5.1. 62 Таблица 5.1 Точное и приближенное решение задачи Коши (5.3) i ti Аналитическое решение y(ti) Метод Эйлера yi Метод Рунге-Кутта второго порядка yi 0 0,5 0,5 0,5 1 1,2 0,72 0,7 0,716667 2 1,4 0,98 0,933333 0,972620 3 1,6 1,28 1,2 1,26788 4 1,8 1,62 1,5 1,60246 5 2 2 1,83333 1,97636 Сравнивая их с результатами аналитического расчета, можно заключить, что приближенные значения содержат существенную погрешность, нарастающую по мере увеличения ti. Для того, чтобы получить этим методом достаточно точное решение, нам потребу- ется значительно уменьшить шаг. Решим теперь эту задачу методом Рунге-Кутта второго порядка. 1. Для первого узла (i = 1) 0 0 0 0,52 2 1 1 yk t     ; *1 0 0 0,5 1 0,2 0,7y y k t       , * * 1 1 1 0,72 2 1,16667 1,2 yk t     ,    *1 0 0 11 10,5 1 1,16667 0,2 0,7166672 2y y k k t         ; 2. Для второго узла (i = 2) 1 1 1 0,7166672 2 1,19445 1,2 yk t     , * 2 1 1 0,716667 1,19445 0,2 0,955557y y k t       , * * 2 2 2 0,9555572 2 1,36508 1,4 yk t     ,    *2 1 1 21 10,716667 1,19445 1,36508 0,2 0,972620.2 2y y k k t         63 Аналогично рассчитываются значения yi в остальных узлах (табл. 5.1). Анализируя полученные результаты, можно заключить, что погрешность решения при использовании метода Рунге-Кутта второго порядка существенно уменьшилась. 5.4. Применение методов Рунге-Кутта к решению систем дифференциальных уравнений первого порядка Пусть имеется задача Коши для системы n дифференциальных уравнений первого порядка  , , ,...y f t y z , y(t0) = y0  , , ,...z g t y z , z(t0) = z0 ……… , (5.4) где t – независимая переменная или аргумент; y = y(t), z = z(t), … – искомые функции аргумента t. Решая систему n дифференциальных уравнений с n неизвестны- ми функциями, мы находимся в пространстве (n + 1) измерений. Положение точки в таком пространстве определяется (n + 1) числа- ми – моментом времени t и значениями функций y и z. Через каж- дую точку пространства проходит n интегральных кривых y = y(t), z = z(t), … . К каждой из них можно провести свою касательную и определить тангенсы углов наклона через правые части соответ- ствующих уравнений: k = f(t, y, z, …), l = g(t, y, z, …), … . Остается только обобщить формулы (5.1) и (5.2). Получаем: 1) для метода Эйлера (Рунге-Кутта первого порядка) yi = yi-1 + ki-1∆t, zi = zi-1 + li-1∆t (i = 1, 2, …), ………, (5.5) где ki-1 = f(ti-1, yi-1, zi-1, …) – тангенс угла наклона касательной к графику y в точке (ti-1, yi-1, zi-1, …); li-1 = g(ti-1, yi-1, zi-1, …) – то же к графику z. 64 2) для метода Рунге-Кутта второго порядка yi = yi-1 + 1/2(ki-1 + ki*)∆t zi = zi-1 + 1/2(li-1 + li*)∆t, (i = 1, 2, …), ……… , (5.6) где ki-1 и li-1 – то же, что и в формуле (5.5); ki* = f(ti, yi*, zi*, …) – тангенс угла наклона касательной к графику y в точке (ti, yi*, zi*, …), причем yi* = yi-1 + ki-1∆t, zi* = zi-1 + li-1∆t; li* = g(ti, yi*, zi*, …) – то же, что и ki*, но к графику z. Формулы численного интегрирования как отдельных дифферен- циальных уравнений, так и их систем не являются точными. На каждом шаге вычислительного процесса они дают некоторую погрешность, способную накапливаться и зависящую от порядка метода, величины шага, свойств самого решения. Поэтому встает вопрос не просто о численном решении задачи Коши, а о получении результата заданной точности. Существует несколько подходов к решению этой проблемы. Самый простой, но и самый трудоемкий подход состоит в следую- щем: пусть требуется получить решение задачи Коши для диффе- ренциального уравнения на отрезке [t0, tk] с m верными значащими цифрами. Разобьем отрезок [t0, tk] на две части, найдем шаг ∆t = (tk – t0)/2 и зафиксируем на отрезке узлы t1 и t2 = tk. Решив зада- чу Коши методом первого или второго порядка, получим значения  0 1y и  02y . После этого уменьшим шаг вдвое и, зафиксировав узлы t1, t2, t3, t4 = tk, снова решим ту же задачу и получим значения  11y ,  1 2y ,  1 3y ,  1 4y . Проверяем условие (1) (0) (1) 4 2 4y y y     . Это значит, что относительное расхождение между реше- ниями на конце отрезка сравнивается с малой величиной, которую можно взять, например, равной 5·10–(m + 1), где m – требуемое коли- чество значащих цифр. Если это условие выполняется, то решение 65  1 1y ,  1 2y ,  1 3y ,  1 4y можно считать искомым. В противном случае шаг ∆t вновь уменьшается вдвое, из решения задачи находятся  21y , …,  28y , проверяется условие (2) (1) (2)8 4 8y y y     и т. д. 5.5. Пример решения задачи Коши для системы двух линейных дифференциальных уравнений первого порядка Пусть требуется решить задачу Коши y yt z  ,y(0) = 1, z t y  ,z(0) = –0,2, (5.7) получив значения yi, zi с четырьмя верными значащими цифрами на отрезке [0, 0,5] с шагом ∆t = 0,1. Для первых двух узлов имеем: Метод Эйлера (Рунге-Кутта первого порядка) 1. Для первого узла (i = 1) k0 = f(t0, y0, z0) = t0y0 + z0 = 0·1 – 0,2 = –0,2, l0 = g(t0, y0, z0) = t0 + y0 = 0 + 1 = 1, y1 = y0 + k0∆t = 1 – 0,2·0,1 = 0,98, z1 = z0 + l0∆t = –0,2 + 1·0,1 = –0,1; 2. Для второго узла (i = 2) k1 = f(t1, y1, z1) = t1y1 + z1 = 0,1·0,98 – 0,1 = –0,002, l1 = g(t1, y1, z1) = t1 + y1 = 0,1 + 0,98 = 1,08, y2 = y1 + k1∆t = 0,98 – 0,002·0,1 = 0,9798, z2 = z1 + l1∆t = –0,1 + 1,08·0,1 = 0,008. Метод Рунге-Кутта второго порядка 1. Для первого узла (i = 1) k0 = f(t0, y0, z0) = t0y0 + z0 = 0·1 – 0,2 = –0,2, l0 = g(t0, y0, z0) = t0 + y0 = 0 + 1 = 1, y1* = y0 + k0∆t = 1 – 0,2·0,1 = 0,98, z1* = z0 + l0∆t = –0,2 + 1·0,1 = –0,1, 66 k1* = f(t1, y1*, z1*) = t1y1* + z1* = 0,1·0,98 – 0,1 = –0,002, l1* = g(t1, y1*, z1*) = t1 + y1* = 0,1 + 0,98 = 1,08, y1 = y0 + 1/2·(k0 + k1*)∆t = 1 + 1/2·(–0,2 – 0,002)·0,1 = 0,9899, z1 = z0 + 1/2·(l0 + l1*)∆t = –0,2 + 1/2·(1 + 1,08)·0,1 = –0,096; 2. Для второго узла (i = 2) k1 = f(t1, y1, z1) = t1y1 + z1 = 0,1·0,9899 – 0,096 = 0,00299, l1 = g(t1, y1, z1) = t1 + y1 = 0,1 + 0,9899 = 1,0899, y2* = y1 + k1∆t = 0,9899 + 0,00299·0,1 = 0,990199, z2* = z1 + l1∆t = –0,096 + 1,0899·0,1 = 0,01299, k2* = f(t2, y2*, z2*) = t2y2* + z2* = 0,2·0,990199 + 0,01299 = 0,21103, l2* = g(t2, y2*, z2*) = t2 + y2* = 0,2 + 0,9990199 = 1,190199, y2 = y1 + 1/2·(k1 + k2*)∆t = 0,9899 + 1/2·(0,00299 + 0,211032)·0,1 = 1,000601, z2 = z1 + 1/2·(l1 + l2*)∆t = –0,096 + 1/2·(1,0899 + 1,190199)·0,1 = 0,018005. Значения yi и zi в остальных узлах рассчитываются аналогично. Полученные результаты показаны в табл. 5.2 и на рис. 5.2. Как вид- но, значения, полученные двумя методами, отличаются между собой и, конечно же, от истинных значений. Таблица 5.2 Точное и приближенное решение задачи Коши (5.7) i ti Метод Эйлера Метод Рунге-Кутта второго порядка yi zi yi zi 0 0 1 –0,2 1 –0,2 1 0,1 0,98 –0,1 0,9899 –0,096 2 0,2 0,9798 0,008 1,000601 0,018005 3 0,3 1,000196 0,12598 1,033747 0,144156 4 0,4 1,0428 0,256 1,091921 0,284802 5 0,5 1,110112 0,40028 1,178801 0,442602 Повысить точность полученного решения можно, уменьшив шаг. Однако это потребует огромной вычислительной работы, особенно Повысить точность полученного решения можно, уменьшив шаг. Однако это потребует огромной вычислительной работы, особенно для метода первого порядка. Для сокращения числа арифметиче- ских действий надо стремиться к использованию метода более 67 высокого порядка. Хотя на каждом шаге он требует выполнения большего числа действий, их общее количество для достижения результата той же точности оказывается меньшим. Еще большего снижения трудоемкости можно добиться при использовании метода четвертого порядка (необходимые формулы можно найти в литера- туре по численным методам). -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0 0,1 0,2 0,3 0,4 0,5 y(Э) z(Э) y(РК) z(РК) t y, z Рис. 5.2. Интегральные кривые искомых функций в задаче Коши (5.7) 5.6. Применение методов Рунге-Кутта к решению дифференциальных уравнений второго и старшего порядков Как следует из вышеизложенного, методы Рунге-Кутта предна- значены для решения отдельных дифференциальных уравнений первого порядка и систем таких уравнений. Если же возникает необходимость в решении уравнений второго или старших поряд- ков, то исходную задачу Коши надо предварительно преобразовать, заменив уравнение n-го порядка с одной неизвестной функцией – системой n уравнений первого порядка с n не- известными функция- ми. Соответственно этому изменяются и начальные условия. В качестве первого примера рассмотрим задачу Коши, которой описываются вынуж- денные колебания си- стемы с одной степе- m y=y(t)` P=P(t) Рис. 5.3. Динамическая система с одной степенью свободы 68 нью свободы. Пусть в произвольной точке однопролетной балки (рис. 5.3) установлен электромотор, деревообрабатывающий станок либо иное механическое устройство с массой m. От веса этого устройства балка изгибается, занимая положение статического рав- новесия, показанное на рис. 5.3 штрих-пунктирной линией. Когда устройство работает, оно создает динамическую нагрузку P = P(t), под действием которой масса (а с ней и балка) совершает колебания вверх-вниз около положения статического равновесия (штриховые линии). Масса балки обычно мала по сравнению с массой установленно- го на ней устройства. Пренебрегая ею, мы получаем систему с одной степенью свободы. Из динамики известно, что вынужденные колебания такой системы описываются следую- щей задачей Коши: ( )my by ry P t    , y(t0) = y0, 0 0( )y t y  , (5.8) где m – масса устройства, кг; b – коэффициент сопротивления, величина которого зависит, главным образом, от того, в какой среде происходят колебания (воздухе, воде, машинном масле и т. д.), кг/с; r – коэффициент жесткости балки в точке расположения устрой- ства; численно равен той силе, которую нужно приложить вместо массы, чтобы балка в той точке получила единичное перемещение, кг/см2; P(t) – динамическая нагрузка – сила, меняющая во времени свою величину по некоторому известному закону H; y0, (0)y – начальное перемещение (м) и начальная скорость (м/с) массы, т. е. перемещение и скорость в момент времени t = t0. Таким образом, имеем задачу Коши для уравнения второго порядка. Ее решение следует начать с преобразования. Для этого вводится обозначение y z , которое одновременно рассматривает- ся как первое уравнение преобразованной задачи. Второе уравнение получается из заданного, проведя в нем замены y z  и y z  mz bz ry P t   или  1 ( )z P t bz ry m    . 69 Первое начальное условие y(t0) = y0 полностью подходит новой задаче, а второе 0 0( )y t y  требует изменения: z(t0) = z0, где z0 чис- ленно равно 0y . Окончательно преобразованная задача выглядит следующим образом: y z , y(t0) = y0,  1 ( )z P t bz ry m    , 0 0 0( )z t z y   . (5.9) Если сравнить ее с записью задачи Коши для задачи (5.4, п.5,4), то в данном случае f(t, y, z) = z,  1( , , )g t y z P t bz ry m       . Теперь ее можно решать методами Рунге-Кутта. Следует особо отметить, что при численном решении этой зада- чи (5.9) мы вынуждены находить не только узловые значения yi, как того требует исходная задача, но и zi. Это как бы побочный продукт, вызванный нашим неумением численно решать уравнения второго и старшего порядков. Второй пример – жестко защемленная консоль постоянного сечения, загруженная произвольной поперечной нагрузкой q(x) (рис. 5.4). q (x ) l x y= y (x ) x Рис. 5.4. К задаче Коши об изгибе консоли 70 Задача состоит из обыкновенного дифференциального уравнения четвертого порядка и четырех начальных условий: IV ( )q xy EJ  , (5.10) y(0) = 0, y'(0) = 0, M(0) = M0, Q(0) = Q0, где y = y(x) – функция прогиба балки; EJ = const – жесткость балки на изгиб; M0, Q0 – изгибающий момент и поперечная сила в начальном се- чении x0 = 0. Чтобы преобразовать уравнение четвертого порядка к системе четырех уравнений первого порядка, введем три вспомогательные функции:  – угол поворота сечения балки, M – изгибающий момент, Q – поперечную силу, а затем воспользуемся дифференци- альными зависимостями, известными из курса сопротивления материалов  = y', M = –EJy" = –EJ', Q = M' (здесь учтено, что положительными являются: прогиб – вниз, угол поворота – по часовой стрелке, изгибающий момент – растягиваю- щий нижние волокна, поперечная сила – вращающая отсеченную часть по часовой стрелке). В итоге получаем следующую преобразованную задачу Коши y' = , y(0) = 0; ' = –M/EJ, (0) = 0; M' = Q, M(0) = M0; ' = –q, Q(0) = Q0. Эту задачу также можно решить рассмотренными в п. 5.4 методами. 71 5.7. Решение задачи Коши в системе Mathematica Численные методы решения дифференциальных уравнений реа- лизуются с помощью функции NDSolve: NDSolve[f,y[x],{x,xmin,xmax}] NDSolve[{f1,f2,...,y1(x0),y2(x0),...},{y1[x],y2[x], ...},{x,xmin,xmax}] f – дифференциальное уравнение и начальные условия; fi – i-е уравнение системы дифференциальных уравнений; y[x] – искомая функция, yi[x] – i-я искомая функция системы диффе- ренциальных уравнений; yi(x0) – i-е начальное условие; xmin, xmax – минимальное и максимальное значения независимой переменной; x – аргумент искомой функции. Рассмотрим решение рассмотренных в лабораторной работе примеров в системе Mathematica. Решение задачи (5.3): 72 Решение задачи (5.7): 5.8. Выполнение индивидуальных заданий Используя трехзначный шифр, выбрать из табл. 5.3 данные к задаче Коши для динамической системы (уравнение 5.8, п. 5.6) с одной степенью свободы. Решить ее методом Рунге-Кутта первого и второго порядка. Результаты представить в виде таблиц и графи- ков. Решить задание, используя систему Mathematica и сравнить полученные результаты. 73 Таблица 5.3 Исходные данные к заданию Шифр 1 цифра шифра 2 цифра шифра 3 цифра шифра № m b r P(t) y0 y0 t0 tk ∆t 1 2 3 4 5 6 7 8 9 10 1 3,0 0,1 2,4 sin(t) + 6,4 0,1 1,0 0 1 0,1 2 3,5 0,3 2,7 sin(2t) – t –0,2 –3,0 0,1 2,1 0,2 3 4,0 0,5 3,0 4sin(2t) 0,3 5,0 0,2 3,2 0,3 4 4,5 0,7 3,3 6cos(3t) –0,4 –7,0 0,3 4,3 0,4 5 5,0 0,9 3,6 5(t 2 – cos(t)) 0,5 1,5 0,4 5,4 0,5 6 5,5 0,2 3,9 6t – cos(2t) –0,6 –2,6 0,5 2,0 0,15 7 6,0 0,4 4,2 7sin(3t) 0,7 8,0 0,6 3,1 0,25 8 6,5 0,6 4,5 8cos(2t) –0,8 –6,0 0,7 4,2 0,35 9 7,0 0,8 4,8 0,5t 2 + 7,3 0,9 4,0 0,8 4,8 0,4 0 7,5 1,0 5,1 –0,2t + 6,5 –1,0 2,0 0,9 5,4 0,45 ЛАБОРАТОРНАЯ РАБОТА № 6 РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ ДЛЯ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ МЕТОДОМ КОНЕЧНЫХ РАЗНОСТЕЙ Цель работы: применить метод конечных разностей к решению краевой задачи для линейных дифференциальных уравнений. Состав работы: 1. Применение метода конечных разностей к решению линейных дифференциальных уравнений любого порядка; 2. Пример решения линейного дифференциального уравнения с заданными краевыми условиями; 3. Формулировка краевых задач для сжато-изогнутого стержня; 4. Решение краевых задач в системе Mathematica; 5. Выполнение индивидуальных заданий. 6.1. Применение метода конечных разностей к решению линейных дифференциальных уравнений любого порядка Краевыми задачами называются задачи, в которых для решения линейных дифференциального уравнения дополнительные условия 74 задаются при двух значениях независимой переменной на концах рассматриваемого участка. Конкретные краевые условия (иногда их называют граничные условия) выбираются из физической поста- новки задачи. Для одного уравнения первого порядка краевую зада- чу поставить нельзя, так как в этом случае для выделения частного решения требуется только одно дополнительное условие, а в крае- вой задаче два конца, и, следовательно, минимальное количество дополнительных условий равно двум т. е., краевую задачу можно сформулировать только для дифференциальных уравнений второго и старшего порядков. Для уравнений второго порядка возможны следующие прос- тейшие варианты краевых условий: 1) на двух концах отрезка зада- ются значения искомой функции; 2) на одном из концов отрезка задается значение искомой функции, на втором – значение первой производной. Основная идея метода конечных разностей (МКР) состоит в том, что он сводит решение краевой задачи к решению системы алгебраических уравнений относительно значений искомой функ- ции (решения) на заданном множестве точек. Это достигается заме- ной производных, входящих в дифференциальное уравнение и в краевые условия, их конечными разностями. Рассмотрим аппроксимацию первой производной для функции y = f(x), заданной таблично в виде y = y1, y2, … , yn-1, yn при x = x1, x2, … , xn-1, xn (рис. 6.1). Считаем, что узлы равноудалены друг от друга. Тогда шаг ∆x (разность между соседними значениями ар- гумента) будет постоянным. x1 x2 xnxn-1 yn y1 yi-1 yi yi+1 a bx x xi xi+1xi-1 Рис. 6.1. Фрагмент отрезка [a, b] с нанесенными на нем узлами и узловыми ординатами искомой функции 75 В зависимости от способа вычисления конечных разностей получают разные формулы для вычисления первой производной (первая разность) в одной и той же точке: левая разность 1i ii y yy x    ; (6.1) правая разность 1i ii y yy x     ; (6.2) центральная разность 1 1 2 i i i y yy x     . (6.3) Конечно-разностные выражения для старших производных получаются по принципу «разность от разности». Можно получить формулы для – второй производной (вторая разность)    1 11 1 1 2 / / 2( ) i i i ii i i i ii i y y x y y xy y y y yy y x x x                   ;(6.4) – третьей производной (третья разность)    2 21 2 2 11 1 2 / 2 / 2 2 i i i i i ii i i y y y x y y y xy yy x x                 2 1 1 2 3 2 2 ; 2 i i i iy y y y x         (6.5) – четвертой производной (четвертая разность) IV 1 1 2 1 1 1 1 2 2 4 2 ( 2 ) 2( 2 ) ( 2 )i i i i i i i i i i i i i y y y y y y y y y y y yy x x                       2 1 1 2 4 4 6 4i i i i iy y y y y x         (6.6) и т. д., являющиеся формулами численного дифференцирования, точность которых существенно зависит от ∆x. 76 Пусть имеется краевая задача, состоящая из линейного диффе- ренциального уравнения второго или старших порядков и краевых условий на концах отрезка [a, b]. Требуется найти значения искомой функции в узлах, равномерно расположенных на отрезке.Решение краевой задачи методом конечных разностей можно разделить на четыре основных этапа. Этап 1. Подготовка к решению задачи. Отрезок [a, b] разбивается на m частей (шагов) длиной ∆x = (b – a) / m. Точки деления (узлы) нумеруются слева направо: например 1, 2, … , m. Узлы 1 и m называются контурными, а узлы 2, 3, … , m–2, m–1 – внутренними. В некоторых случаях приходится вводить еще и законтурные узлы, расположенные за пределами от- резка. Неизвестными задачи выступают узловые значения искомой функции во внутренних узлах. Если значение искомой функции в одном из контурных узлов не задано краевыми условиями, то в число неизвестных надо включить узловое значение функции в соответствующем контурном узле. Этап 2. Составление прямоугольной СЛАУ. В решаемом дифференциальном уравнении аргумент x заменяет- ся узловым значением xi, искомая функция y(x) – узловым значением yi = y(xi), производные искомой функции ( )y x , ( )y x , … – соответ- ствующими конечными разностями iy , iy , … .Если в уравнении содержится некоторая неизвестная функция аргумента x, например, f(x), то она также заменяется узловым значением fi = f(xi). В результате таких замен дифференциальное уравнение превра- щается в алгебраическое, называемое конечно-разностным для про- извольного узла i или шаблоном. По этому шаблону составляются уравнения :во-первых, во всех внутренних узлах (i = 2, … , m–1), во-вторых, для того контурного узла 1 или m, в котором не задана величина ya или yb . В совокупности составленные уравнения образуют СЛАУ, име- ющую прямоугольную структуру, т.к. число неизвестных в ней пре- вышает количество уравнений. Причем превышает на величину, равную порядку решаемого уравнения или, что одно и то же, коли- честву краевых условий. Этап 3. Приведение прямоугольной СЛАУ к квадратному виду. Такое приведение может быть выполнено, по крайней мере, 77 двумя способами. Первый состоит в том, что к составленным добавляют недостающие уравнения, в качестве которых использу- ются краевые условия. Второй способ заключается в преобразова- нии составленных уравнений с помощью краевых условий, в ре- зультате чего число неизвестных сокращается до числа уравнений. Первый способ более простой как при ручном, так и автоматизиро- ванном формировании системы, но ее порядок оказывается выше, чем при использовании второго способа. Этап 4. Решение квадратной СЛАУ. Используя любой метод (например, метод Гаусса, метод прогон- ки и т. д.) систему решают и получают приближенные узловые зна- чения искомой функции. Приближенность результатов решения краевой задачи методом конечных разностей обусловлена влиянием ошибки метода и ошибки округления. Первая возникает вследствие замены производных приближенными конечно-разностными выра- жениями, а вторая – главным образом при решении СЛАУ. Сделать первую ошибку малой можно, дробя шаг (при этом возрастают размеры СЛАУ), а вторую – переходя к машинной арифметике удвоенной точности. 6.2. Пример решения линейного дифференциального уравнения с заданными краевыми условиями методом конечных разностей Настоящий пример не ставит своей целью получение результа- тов высокой точности, а служит для понимания МКР. Поэтому шаг ∆x в задаче взят довольно крупным. Задача. Решить методом конечных разностей краевую задачу (число шагов m = 4) x2y" – 2xy' + 2y = x3 ; y(1) = 1,1; y(2) = –0,5. Решение. 1. Отрезок [a, b], где a = 1; b = 2, разбиваем на m = 4 частей (шагов) длиной ∆x = (2 – 1)/4 = 0,25. 78 Узлы нумеруем слева направо от 1 до 5. Неизвестными являются величины y2, y3, y4 (рис. 6.2). Известные величины – y1 = ya = 1,1 (при x1 = 1) и y5 = yb = –0,5 (при x5 = 2). x1,25 1,5 1,75 (x1,y1) x=0,25 Интегральная кривая, отвечающая точному решению 0,5 1,0 1,0 2,0 (x2,y2) (x3,y3) (x4,y4) (x5,y5) Рис. 6.2. К решению краевой задачи 2. Используя конечно-разностные выражения (6.3) и (6.4), преобразуем исходное дифференциальное уравнение 2 31 1 1 1 2 2 2 2 2 i i i i i i i i i y y y y yx x y x xx         . Умножая все уравнения на ∆x2 и приводя подобные слагаемые, получим 2 2 2 2 3 2 1 1( ) 2( ) ( )i i i i i i i i ix x x y x x y x x x y x x           . Учитывая, что ∆x = 0,25; x2 = 1,25; x3 = 1,5; x4 = 1,75, составляем систему трех алгебраических уравнений i = 2: (1,252 + 1,25·0,25)y1 + 2(0,252 – 1,252)y2 + (1,252 –1,25·0,25)y3 = 1,253·0,252; i = 3: (1,52 + 1,5·0,25)y2 + 2(0,252 – 1,52)y3 + (1,52 –1,5·0,25)y4 = 1,53·0,252; i = 4: (1,752 + 1,75·0,25)y3 + 2(0,252 – 1,752)y4 + (1,752 –1,75·0,25)y5 = 1,753·0,252. 79 Подсчитывая значения коэффициентов, записываем СЛАУ. 1 2 3 2 3 4 3 4 5 1,875 3 0,9375 0,122070, 2,625 4,375 1,875 0, 210938, 3,5 6 2, 625 0,334961. y y y y y y y y y          Получили систему трех уравнений с пятью узловыми значения- ми искомой функции: y1, y2, y3, y4, y5. 3. Для того, чтобы избавится от «лишних» переменных, восполь- зуемся вторым способом, описанном в этапе 3 п. 6.1. Вместо y1 в первом уравнении подставляем значение ya = 1,1, а вместо y5 в третьем уравнении – значение yb = –0,5. Перенося числа 1,875·1,1 и 2,625·(–0,5) в правую часть соответ- ственно первого и третьего уравнений и алгебраически складывая их со стоящими там свободными членами, окончательно получим квадратную СЛАУ третьего порядка и ленточной (трехдиагональ- ной) структуры: 2 3 2 3 4 3 4 3 0,9375 1,94043, 2,625 4,375 1,875 0, 210938, 3,5 6 1, 64746. y y y y y y y          (6.7) 4. Решая полученную систему линейных уравнений любым из известных методов, находим y2 = 0,770252, y3 = 0,395015, y4 = –0,0441510. Тогда исходная задача имеет следующее решение y = (1,1; 0,770252; 0,395015; –0,0441510; –0,5) при x = (1; 1,25; 1,5; 1,75; 2). По этим значениям на рис. 6.2 построена сплошная линия. 80 Рассмотренная краевая задача имеет точное аналитическое решение y(x) = 0,5x(6,9 – 5,7x + x2), которому отвечает пунктирная линия на рисунке 6.2. При этом y(x2) = 0,835938, y(x3) = 0,45, y(x4) = –0,0109375. Как и ожидалось, метод конечных разностей при столь большом шаге дает только качественно правильное решение. 6.3. Формулировка краевых задач для сжато-изогнутого стержня Рассмотрим упругий стержень длиной l постоянного поперечно- го сечения с жесткостью на изгиб EJ = const, загруженный произ- вольной поперечной нагрузкой q(x) и сжатый силой P, не превы- шающей своего критического значения Pcr.. Как известно из курса сопротивления материалов, работа такого стержня описывается обыкновенным дифференциальным уравнением четвертого порядка EJyIV + Py" = q(x), (6.8) где y = y(x) – функция прогибов оси стержня, и четырьмя краевыми условиями, отражающими характер опирания стержня по двум его концам. Приведем формулировку краевых условий и формулы критиче- ских сил для шести схем стержня. Во всех схемах начало координат совмещено с левым концом стержня. Ось x направлена слева напра- во вдоль оси стержня, так что a = 0, b = l. Ось прогибов y направле- на вниз, такое же направление имеет и нагрузка q(x). 81 Схема 1 P q(x) x y y(x) l y(0) = 0, y(l) = 0, y"(0) = 0,y"(l) = 0 (это значит, что на шар- нирной опоре равны ну- лю прогиб и изгибающий момент (см. также приме- чание после схемы 6)). 2 2cr EJP l   . Схема 2 P q(x) x y y(x) l y(0) = 0,y(l) = 0, y'(0) = 0, y"(l) = 0 (в защемлении равны ну- лю прогиб и угол поворо- та). 2 2(0,7 )cr EJP l   . Схема 3 P q(x) x y y(x) l y(0) = 0,y(l) = 0, y'(0) = 0, y'(l) = 0 (скользящее защемление правого конца стержня при малых прогибах опи- сывается так же, как и обычная заделка). 2 2(0,5 )cr EJP l   . 82 Схема 4 P q(x) x y y(x) l y(0) = 0, y'(l) = 0, y"(0) = 0,y'''(l) = 0 (в скользящем защем- лении на правом конце равны нулю угол поворота и поперечная сила). 2 2(2 )cr EJP l   . Схема 5 P q(x) x y y(x) l y(0) = 0, y'(l) = 0, y'(0) = 0, y'''(l) = 0, 2 2cr EJP l   . Схема 6 q(x) l x y y(x) P y(0) = 0, y"(l) = 0, y'(0) = 0, –EJy'''(l) = Py'(l), 2 2(2 )cr EJP l   . 83 Примечание. При формулировке краевых условий учтены дифференциальные зависимости при изгибе, известные из курса сопротивления материалов: y' = , M = –EJy", Q = –EJy'''. Отсюда равенство нулю угла поворота  в заделке x = 0, x = l записывается как y'(0) = 0 или y'(l) = 0. Равенство нулю изгибающе- го момента в шарнире и на концах x = 0, x = l отражается в форму- лировке y"(0) = 0 или y"(l) = 0. Аналогично записывается четвертое краевое условие в схе- мах 4 и 5: y'''(l) = 0. Оно отражает факт равенства нулю поперечной силы Q на пра- вом конце балки x = l. Несколько сложнее формулируется краевое условие в схеме 6. При изгибе стержня его сечение поворачивается на угол  (рис. 6.3), а сжимающая сила P сохраняет свое прежнее (горизон- тальное) положение. Разложив ее на два направления (нормальное и касательное к оси деформированного стержня), получаем, что на конце стержня N = –P cos ≈ –P, Q = P sin ≈ P y'. Замены cos ≈ 1, sin ≈ 0 возможны только при малых прогибах. P   Psin Pcos Рис. 6.3. К формулировке четвертого краевого условия в схеме 6 84 6.4. Решение краевых задач в системе Mathematica Для решения линейного дифференциального уравнения с задан- ными краевыми условиями из п. 6.2 воспользуемся функцией системы Mathematica NDSolve (см. ее описание в лабораторной работе № 5) 6.5. Выполнение индивидуальных заданий Задана краевая задача, состоящая из обыкновенного дифферен- циального уравнения второго порядка y" + ky = f(x) и двух краевых условий y(0) = y(l) = 0. Значения k = k*(π/l)2, l и выражение f(x) принимаются по табл. 6.3. Используя метод 85 конечных разностей при m = 4, составить таблицу с результатами и построить график функции y = y(x). При решении системы линей- ных алгебраических уравнений использовать метод прогонки. Решить задание, используя систему Mathematica и сравнить полу- ченные результаты. Таблица 6.3 Исходные данные к заданию Первая цифра шифра k* Вторая цифра шифра l Третья цифра шифра f(x) 1 2 3 4 5 6 0 0,65 0 7,6 0  20,01 x lxl  1 0,45 1 5,6 1 2 0,15 sin xx ll   2 0,5 2 8,4 2 3 8( 3 ) ( 1) l x l x   3 0,55 3 6,4 3 3 2 225 x l x l  4 0,8 4 6,8 4 2 2ln 1 5 x x ll      5 0,6 5 7,2 5 3 1,5 2x ll   6 0,35 6 6,0 6 31,5 sin x x ll  7 0,75 7 8,8 7 2 3 3 2 2 100 lx x l l   8 0,7 8 8,0 8     2 3 15ln 2 2 x l x l   9 0,4 9 5,2 9 3 3 4 3x ll     86 Указания к работе. 1. Работа шарнирно опертого по обоим концам стержня длиной l постоянной жесткости на изгиб EJ, загруженного произвольной по- перечной нагрузкой и сжатой силой P (рис. 6.4), описывается обык- новенным дифференциальным уравнением ( )P M xy y EJ EJ    и краевыми условиями y(0) = y(l) = 0, где y = y(x) – функция прогибов стержня; M(x) – изгибающий момент в произвольном сечении x. Введя обозначения k = P/EJ и f(x) = M(x)/EJ, получаем заданное уравнение. 2. Разделив пролет стержня на 4 части (m = 4), надо отметить на оси x пять узлов 1, 2, …, 5, отстоящих друг от друга на шаг ∆x = l/4. При составлении конечно-разностного шаблона решаемого уравне- ния правая часть заменяется узловым значением fi = f(xi), последова- тельно полагая i равным 2, 3, 4, получают три уравнения с пятью неизвестными: y1, y2, y3, y4, y5. Чтобы сделать из этих уравнений квадратную СЛАУ, надо положить в них (в соответствии с краевы- ми условиями) y1 = y5 = 0. Решив эту систему, находят приближен- ные значения прогибов стержня. Выбрав масштаб по своему усмот- рению, нужно построить график прогибов стержня (положительные значения прогиба откладываются вниз от оси). P q(x) x y y(x) l x Рис. 6.4. Расчетная схема сжато-изогнутого стержня 87 3. В лабораторной работе № 2 для решения систем линейных алгебраических уравнений использовался метод Гаусса. Вместе с тем, для решения СЛАУ специального вида, например, с трехдиа- гональной структурой (6.7), можно использовать метод прогонки. Рассмотрим СЛАУ 3 (n = 3) порядка следующего вида: 1 1 1 2 1 2 1 2 2 2 3 2 3 2 3 3 3 , , b x c x d a x b x c x d a x b x d        Метод состоит из двух ходов и напоминает метод Гаусса. Прямой ход заключается в поиске коэффициентов P и Q по следующим формулам: P1 = 0 , Q1 = 0, 1 ii i i i cP a P b   , 1 i i i i i i i d a QQ a P b   , 1,i n . Во время обратного хода непосредственно определяются значе- ния неизвестных в обратной последовательности: xn = Qn + 1, 1 1 1i i i ix Q P x    , 1, 2, 1i n n    . 88 ЛАБОРАТОРНАЯ РАБОТА № 7 ЧИСЛЕННЫЕ МЕТОДЫ ОПТИМИЗАЦИИ. ГРАФИЧЕСКИЙ МЕТОД РЕШЕНИЯ ЗАДАЧ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ Цель работы: изучить графический метод решения задач ли- нейного программирования. Состав работы: 1. Общие сведения о задачах оптимизации; 2. Графическое решение ЗЛП; 3. Решение задач линейного программирования в системе Mathematica; 4. Выполнение индивидуальных заданий. 7.1. Общие сведения о задачах оптимизации В настоящей работе будем рассматривать задачи оптимизации, имеющие следующую структуру: 1) задана функция цели Z = f (x1, x2, … xn), которая принимает максимальное или минимальное значение; 2) присутствуют m ограничений 1 2( , , , )i n ig x x x b        , i = 1, … m; 3) задано условия неотрицательности переменных x1 ≥ 0, x2 ≥ 0, … , xn ≥ 0. Задача оптимизации в такой постановке называется условной. Вторая и третьи части могут отсутствовать, тогда задача называется безусловной. Решить задачу оптимизации – это значит найти такие значе- ния переменных x1, x2, … , xn, которые удовлетворяют ограничени- ям, условиям неотрицательности и сообщают функции цели макси- мальное или минимальное значение. 89 Совокупность указанных значений переменных x1, x2, … , xn будем называть оптимальным планом. Задачи оптимизации описывают явления и процессы в различ- ных областях науки, техники, экономики. Большую роль они игра- ют и в строительстве (оптимальное проектирование конструкций, транспортная задача и другие.) Пусть требуется найти оптимальные размеры b, h прямоугольно- го сечения простой балки пролетом l, загруженной равномерно распределенной нагрузкой интенсивности q (рис. 7.1). Модуль упругости материала E, расчетное сопротивление на растяжение- сжатие R. Допустимый прогиб балки равен c. В качестве функции цели возьмем объем материала, а в качестве ограничений – условия прочности и жесткости. Рис. 7.1. К формулировке задачи оптимизации однопролетной балки Объем материала балки V =bhl должен быть минимальным. Но поскольку пролет l – известная (заданная) величина, то в каче- стве функции цели можно принять площадь сечения. И тогда, вводя обозначения V = Z, b = x1, h = x2, получаем Z = x1x2 → min. Наибольшие нормальные напряжения возникают в фибровых волокнах на середине балки, где действует наибольший изгибаю- щий момент. Причем эти напряжения maxmax M W   90 не должны превышать расчетного сопротивления R. Учитывая, что 2max / 8M ql , а момент сопротивления сечения W = bh2/6, получаем 2 2 6 8 ql R bh   , что после преобразований и замены b = x1, h = x2 дает первое ограничение 2 2 1 2 3 4 qlx x R  . Наибольший прогиб в балке возникает также в середине пролета 45 384 qlf EJ  . Он не должен превзойти допускаемого прогиба c, т. е. 45 384 ql c EJ  . Момент инерции прямоугольного сечения J = bh3/12. С учетом этого выражения второе ограничение примет вид 3 1 2 5 32 qlx x Ec  . Формулируя задачу оптимизации, нужно также учесть, что ис- комые величины b и h не должны принимать отрицательных значе- ний. В итоге получаем задачу оптимизации Z = x1x2 → min, 91 2 2 1 2 3 4 qlx x R  , 3 1 2 5 32 qlx x Ec  , x1 ≥ 0, x2 ≥ 0. Раздел математики, изучающий задачи оптимизации и разраба- тывающий методы их решения, называется математическим программированием (иногда в литературе встречается термин «математическое планирование»). В зависимости от вида функции цели и ограничений математи- ческое программирование делится на линейное, нелинейное, цело- численное, стохастическое, динамическое и т. д. В задачах линейного программирования (ЗЛП) функции цели и ограничения линейны относительно x1, x2, … xn: Z = c1x1 + c2x2+ … +cnxn → max (или min), a11x1 + a12x2 + … +a1nxn ≥ b1, a21x1 + a22x2 + … +a2nxn ≥ b2, ……… am1x1 + am2x2 + … +amnxn ≥ bm, x1 ≥ 0, x2 ≥ 0, … , xn ≥ 0. В рассмотренном выше примере была сформулирована задача нелинейного программирования (ЗНП) – в этих задач функция цели или хотя бы одно из ограничений нелинейны относительно переменных. Эти задачи настолько сложны, что до сих пор не уда- лось разработать общие методы (подобно симплекс-методу в линейном планировании), которые гарантировали получение оп- тимального решения любой ЗНП. Имеются методы решения лишь отдельных специальных классов задач и, прежде всего, задач с выпуклой функцией цели и выпуклыми ограничениями. В задачах целочисленного программирования неизвестные мо- гут принимать только целочисленные значения. Если в целевой функции или в функциях, определяющих область возможных изме- 92 нений переменных, содержатся случайные величины, то такая зада- ча относится к задаче стохастического программирования. Зада- ча, процесс нахождения решения которой является многоэтапным, относится к задаче динамического программирования. В настоящем практикуме рассматриваются ЗЛП и два метода их решения – графический и симплекс-метод. 7.2. Графическое решение ЗЛП Решение таким методом выполняется для ЗЛП с двумя неиз- вестными. Рассмотрим алгоритм графического решения на конкретном примере Z = x1 – 2x2 → min max ; целевая функция 3x1 + 4x2 ≤ 12; 2x1 – x2 ≥ –2; линейные ограничения x1 – 2x2 ≤ 6; x1 ≥ 0, условие неотрицательности где условие неотрицательности наложено только на первую неиз- вестную. По существу, здесь две задачи: одна – это задача миними- зации Z, а вторая – максимизации Z. Решение состоит из четырех этапов. Этап 1. Построение области допустимых решений. Областью допустимых решений (ОДР) в задаче с двумя неиз- вестными называется пересечение полуплоскостей, определяемых ограничениями и условиями неотрицательности. Иначе говоря, это область на плоскости x10x2 , состоящая из точек, координаты которых удовлетворяют как всем ограничениям, так и условиям не- отрицательности. Первое ограничение 3x1 + 4x2 ≤ 12 определяет первую полуплос- кость. Чтобы найти ее, нужно сначала построить граничную пря- мую 3x1 + 4x2 = 12. В данном случае она проходит через точки (0,3) и (4,0) (см. рис. 7.2). Влево и вниз от нее идет одна бесконечная полуплоскость, а вправо и вверх – вторая полуплоскость. Какая 93 из них отвечает ограничению 3x1 + 4x2 ≤ 12 ? Для ответа на этот вопрос возьмем какую-нибудь пробную точку. Если граничная пря- мая не проходит через начало координат, то в качестве пробной точки проще всего использовать именно начало координат, т. е. точку (0,0). Подставляя эти координаты в левую часть первого ограничения, получаем 0 < 12. Таким образом, искомая полуплос- кость включает начало координат. Направление полуплоскости по- казываем стрелками, рядом с которыми пишем номер ограничения. 1 2 4 2 4 1 3 3 2 3 4 5 6 7-1-2-3-4 1 3 4 5 -1 -2 -3 -4 2 c c Z= 0 Область допустимых значений (ОДР) A B D 1 x1 x2 Рис. 7.2. Графическое решение ЗЛП с двумя неизвестными ограничения, получаем 0 < 12. Таким образом, искомая полуплос- кость включает начало координат. Направление полуплоскости по- казываем стрелками, рядом с которыми пишем номер ограничения. Аналогично строятся полуплоскости 2 и 3. Условие неотрица- тельности x1 ≥ 0 дает четвертую полуплоскость. Фигура, очерченная четырьмя граничными прямыми, и есть ОДР. Она представляет собой замкнутую выпуклую область. 94 Этап 2. Построение линии уровня Z = 0. Линией уровня некоторой функции называется такая линия, во всех точках которой значение этой функции одинаково. Линия уровня функции цели при Z = 0 – это прямая c1x1 + c2x2 = 0 или, в нашем случае, x1 – 2x2 = 0. Ее график проходит через начало координат (0, 0) и, например, точку с координатами (2, 1). Во всех точках этой линии Z = 0. Этап 3. Построение вектора-градиента (или вектора-анти- градиента). Вектором-градиентом некоторой функции нескольких пере- менных называется вектор, указывающий направление наискорей- шего возрастания функции. Компоненты вектора-градиента равны частным производным от функции по всем ее переменным. В ЗЛП с двумя переменными Z = c1x1 + c2x2 и 1 1 2 2 / / Z x c c Z x c               . В рассматриваемой задаче  1, 2 Tc   . Соответственно вектор-антиградиент указывает направление наискорейшего убывания функции цели 1 1 2 2 / / Z x c c Z x c                  . В нашей задачи  1, 2 Tc   . Чтобы построить c , от начала координат откладываем вдоль x1 первую координату (равную +1), а вдоль x2 – вторую (равную –2). На этих отрезках строим прямоугольник, диагональ которого, снабженная направлением, и есть искомый вектор-градиент. Он перпендикулярен линии уровня функции цели. Аналогично строится – c , который направлен в сторону, противоположную c . Этап 4А. Решение задачи минимизации Z . Если линию уровня перемещать параллельно самой себе в направлении вектора-антиградиента, то функция цели будет убы- вать, принимая отрицательные значения все большего модуля. 95 И это будет продолжаться до тех пор, пока линия уровня не коснет- ся крайней точки ОДР, а именно, вершины А. Координаты этой вершины и являются искомым решением задачи минимизации Z. Для их определения надо решить систему двух уравнений, соответ- ствующих первой и второй граничным прямым 1 2 1 2 3 4 12 2 2, x x x x      откуда x1 = 4/11 = 0,363636, x2 = 30/11 = 2,72727 (эти значения назы- ваются оптимальным планом при решении задачи минимизации Z). Функция цели принимает в оптимальной вершине А значение min Z = 4/11 – 2·30/11 = –56/11 = –5,09091. Этап 4Б. Решение задачи максимизации Z. Чтобы найти оптимальный план в задаче максимизации Z, линию уровня надо перемещать в направлении вектора-градиента. При этом функция цели будет возрастать. Перемещение линии уровня надо проводить до того момента, пока она не коснется край- ней точки ОДР. Но данная задача имеет особенность: линия уровня параллельна одной из сторон ОДР, а именно стороне BD. Поэтому при поступа- тельном движении линия уровня коснется границы ОДР в беско- нечно большем количестве точек – в точках B, D и во всех проме- жуточных точках. Отсюда следует, что задача максимизации Z в данном случае имеет бесконечное множество решений: оптималь- ным планом являются координаты любой точки отрезка BD. Возьмем для определенности в качестве оптимальной вершины точку B. Ее координаты (оптимальный план) видны на чертеже x1 = 0, x2 = –3, max Z = 0 – 2·(–3) = 6. При решении некоторой ЗЛП может оказаться, что одни ограни- чения противоречат другим (рис. 7.3). Так что найти их пересечение не представляется возможным. В этом случае говорят, что ОДР пустая. Такая задача решения не имеет. Встречаются случаи, когда в направлении вектора-градиента или антиградиента ОДР не замкнута (рис. 7.4). Двигаясь в этом направ- лении, достичь крайней точки ОДР невозможно. Такая ЗЛП также не имеет решения. 96 0 x1 x2 0 Z= 0 ОДР c x1 x2 Рис. 7.3. ОДР пустая (ограничения противоречивы) Рис. 7.4. Задача максимизации Z решения не имеет 7.3. Решение задач линейного программирования в системе Mathematica Для поиска максимума и минимума аналитически заданной функции нескольких аргументов при наличии ограничений с помо- щью системы Mathematica можно использовать следующие функции: NMaximize[{f,cons},{x,y,...}]ищет максимум целе- вой функции f в области, определяемой списком ограничений (неравенствами) cons, {x,y,...} – список переменных; Nminimize [{f,cons},{x,y,...}]ищет минимум це- левой функции f в области, определяемой списком ограниче- ний (неравенствами) cons, {x,y,...} – список переменных. Решим рассматриваемый в лабораторной работе пример с помо- щью системы Mathematica 97 Результаты, полученные в системе Mathematica и с помощью графического метода, совпадают. 7.4. Выполнение индивидуальных заданий Решить графически задачу линейного программирования с дву- мя неизвестными. Исходные данные выбрать по табл. 7.1 в соответ- ствии с шифром. Проверить полученные результаты, используя систему Mathematica. Указания к работе. Функция цели в задаче может достигать как максимума, так и минимума. Таким образом, надо решать по суще- ству две ЗЛП. Условия неотрицательности на переменные x1 и x2 не накладываются. Таблица 7.1 Исходные данные к заданию Первая цифра шифра Коэффициен- ты функции цели Вторая цифра шифра Ограничения Третья цифра шифра d 1 2 3 4 5 6 0 с1 = d c2 = 8 0 3x1 + 5x2 ≥ 2(2 + d) 2x1 – x2 ≥ 2 x1 – x2 ≤ 3 0 2,0 1 с1 = –2 c2 = 2d 1 –2x1 + 1,5x2 ≥ 6 x1 ≤ –1 x1 + x2 ≥ 2(2 – d) 3x1 + 2x2 ≤ 12 1 2,5 2 с1 = 1,5d c2 = –3 2 4x1 + 3x2 ≤ –12 2dx2 ≤ 6 x1 ≥ –7 –2x1 + x2 ≤ 4 2 3,0 3 с1 = –4 c2 = –2,5d 3 x1 – 4x2 ≥ 0 –3x1 + 2x2 ≤ 12 x1 ≤ d – 2 3 3,5 4 с1 = 2d c2 = 5 4 x1 ≥ –1 2x1 – 3x2 ≥ 6 x1 + x2 ≥ –5 4x1 + x2 ≤ 20 – d 4 4,0 98 Окончание табл. 7.1 1 2 3 4 5 6 5 с1 = –6 c2 = 2d 5 –x1 – 2x2 ≤ 8 x1 – 2x2 ≥ –6 dx2 ≤ 10 7x1 – 4x2 ≤ 28 5 3,5 6 с1 = d c2 = –7 6 –dx1 + x2 ≤ 0 x1 – 4x2 ≥ 0 x1 – 2x2 ≤ 8 6 3,0 7 с1 = –8 c2 = –3d 7 2x1 + x2 ≥ 6 x1 – x2 ≤ 2 x1 ≥ 3 – d 7 2,5 8 c1 = 6 c2 = –d 8 2x1 – 3x2 ≤ 6 x1 + x2 ≥ 3 –x1 + x2 ≤ 2 x2 ≤ d + 3 8 2,0 9 c1 = –5 c2 = 2d 9 –x1 + 2x2 ≤ 10 x1 + x2 ≥ d 2x1 – x2 ≤ 4 9 4,0 99 ЛАБОРАТОРНАЯ РАБОТА № 8 СИМПЛЕКС-МЕТОД РЕШЕНИЯ ЗАДАЧ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ Цель работы: изучить симплекс-метод, предназначенный для решения задач линейного программирования. Состав работы: 1. Выводы по графическому решению ЗЛП. Идея и основные этапы симплекс-метода; 2. Жордановы исключения; 3. Приведение ЗЛП к каноническому виду; 4. Пример решения ЗЛП симплекс-методом; 5. Выполнение индивидуальных заданий. 8.1. Выводы по графическому решению ЗЛП. Идея и основные этапы симплекс-метода Графическое решение задач с двумя переменными позволяет сделать два важных вывода: 1) ЗЛП может: а) иметь единственное решение; б) иметь бесконечное множество решений; в) не иметь решения. 2) Если ЗЛП имеет решение, то оно достигается в одной или нескольких вершинах ОДР. Распространяя эти выводы на ЗЛП с n переменными, можно предложить, казалось бы, следующий простой алгоритм: надо найти координаты всех вершин ОДР, вычислить во всех вершинах значе- ния функции цели и, сравнив их между собой, выбрать ту из них, где функция цели наибольшая или наименьшая. Но дело в том, что в задаче с n переменными ОДР представляет собой сложный гео- метрический объект, ограниченный уже не отрезками прямых, а так называемыми гиперплоскостями. Число вершин в таком объекте и сложность определения их координат столь велики, что это не под силу сделать даже компьютеру. Из стремления обойти эти трудности и возник универсальный метод решения ЗЛП с любым числом переменных, а именно, 100 симплекс-метод. Суть его заключается в следующем. Вначале, используя специальный алгоритм, находят любую вершину ОДР (ее координаты называются начальным планом). Применяя некото- рые формальные признаки, эту вершину исследуют на оптималь- ность. Если оказывается, что оптимальной она не является, то опять-таки по строго определенным правилам переходят в сосед- нюю вершину, функция цели в которой больше или, по крайней мере, равна тому значению, которое она имела в начальной вер- шине. Новую вершину также проверяют на оптимальность и т. д. Процесс продолжается до тех пор, пока не будет достигнута опти- мальная вершина либо не будет доказано, что поставленная задача решения не имеет. Опыт показывает, что происходит в итоге срав- нительно небольшого числа шагов. Симплекс-метод состоит из трех основных этапов: 1) приведение ЗЛП к каноническому виду и запись ее в началь- ную жорданову таблицу; 2) нахождение начального плана; 3) нахождение оптимального плана. Каждый из этих этапов будет подробно рассмотрен ниже. Однако вначале познакомимся с аппаратом жордановых исключе- ний, на котором базируются второй и третий этапы симплекс- метода. 8.2. Жордановы исключения Пусть имеется система двух уравнений y1 = 2x1 – 5x2 + 4x3, y2 = 8x1 + 2x2 – 3x3. Переменные x1, x2, x3, стоящие в правой части уравнений, назо- вем свободными или независимыми. Они могут принимать любые значения. Переменные y1, y2, находящиеся в левой части, наоборот, будут зависимыми или несвободными. Их значения определяются выбором свободных переменных. Например, если принять x1 = x2 = x3 = 0, то и y1, y2 будут равны нулю. А если x1 = 1, x2 = 2, x3 = 3, то y1 = 4, y2 = 3. 101 Поставим задачу поменять местами какие-нибудь свободные и несвободные переменные, например, x2 и y2. Такая замена и составляет суть жордановых исключений. Разрешим второе уравнение относительно x2 x2 = –4x1 + 0,5y2 + 1,5x3 и подставим это выражение в первое уравнение y1 = 2x1 – 5(–4x1 + 0,5y2 + 1,5x3) + 4x3 = 22x1 – 2,5y2 – 3,5x3. Получили преобразованную систему y1 = 22x1 – 2,5y2 – 3,5x3, x2 = –4x1 + 0,5y2 + 1,5x3. в которой свободными переменными выступают уже x1, y2, x3, а несвободными – y1 и x2. Запишем исходную и преобразованную системы в жордановы таблицы: –x1 –x2 –x3 –x1 –y2 –x3 1 y1 = –2 5 –4 y1 = –22 2,5 3,5  y2 = –8 –2 3 x2 = 4 –0,5 –1,5 1  3 1 2 3 Обе таблицы имеют верх, в котором помещены свободные пере- менные, взятые со знаком минус. Такая форма их представления удобна для последующего использования в симплекс-методе. Обе таблицы имеют левый столбец, в котором располагаются зави- симые переменные. Коэффициенты систем помещены в две строки и три столбца; их номера расположены слева и снизу от таблиц. Соответственно этим номерам будем считать, что коэффициентами исходной таблицы являются a11 = –2, a12 = 5, … ,а преобразован- ной – a*11 = –22, a*12 = –2,5, … . 102 Тот столбец  из которого переносится независимая перемен- ная x2, будем называть разрешающим столбцом. Ту строку  из которого переносится независимая переменная y2, назовем разрешающей строкой. (Номера этой строки и столбца обведены кружком). На пересечении разрешающей строки и разрешающего столбца находится разрешающий элемент a22 = –2 (он выделен). Чтобы проверить, что таблицы составлены правильно, можно от них перейти назад к уравнениям. Для этого коэффициенты любой строки надо умножить на стоящие над ними свободные переменные (с учетом их знака!), произведения сложить, а сумму приравнять соответствующей зависимой переменной. Например, из первой строки исходной таблицы получаем первое уравнение исходной системы y1 = –2·(–x1) + 5·(–x2) – 4·(–x3) = 2x1 – 5x2 + 4x3. Сопоставляя исходную и преобразованную таблицы, можно непосредственно убедиться в справедливости следующих четырех правил, по которым происходит преобразование коэффициентов таблиц при жордановых исключениях 1) разрешающий элемент заменяется обратной величиной a22 = –2; a*22 = 1/a22 = 1/(–2) = –0,5; 2) остальные элементы разрешающей строки делятся на разре- шающий элемент a21 = –8, a*21 = –8/(–2) = 4, a23 = 3, a*21 = 3/(–2) = –1,5; 3) остальные элементы разрешающего столбца делятся на раз- решающий элемент и меняют знак на противоположный a12 = 5, a*12 = –5/(–2) = 2,5; 4) все оставшиеся элементы преобразуются по правилу прямо- угольника (рис. 8.1): 103 а) б) (+) a11 (-) a12 a21 a22 (-) a12 (+) a13 a22 a23 Рис. 8.1. Пояснение правила прямоугольников а) элемент a11 = –2 входит в прямоугольник (рис. 8.1 а) в котором он образует положительную диагональ с разрешающим элементом a22. Вторая диагональ a21 – a12 считается отрицательной. Поэтому * 11 22 12 21 11 22 2 ( 2) 5 ( 8) 22; 2 a a a aa a           б) элемент a13 = –4 входит в другой прямоугольник (рис. 8.1 б), образуя положительную диагональ другого направления, но опять- таки с разрешающим элементом a22. Поэтому * 13 22 12 23 13 22 4 ( 2) 5 3 3,5. 2 a a a aa a         В заключении отметим еще одно важное свойство изложенного варианта жордановых исключений. Независимая переменная, нахо- дясь на верху таблицы, имеет знак минус, а при переходе в левый столбец теряет его. Зависимая переменная, переходя из левого столбца наверх, знак минус приобретает. 8.3. Приведение ЗЛП к каноническому виду Для того, чтобы придать симплекс-методу универсальный харак- тер, необходимо привести решаемую задачу к каноническому виду. ЗЛП считается записанной в каноническом виде, если она удовлетворяет четырем требованиям: 1) Функция цели максимизируется. 2) Ограничения – строгие равенства. 3) Свободные члены ограничений неотрицательны. 4) Условия неотрицательности наложены на все переменные. 104 8.4. Пример решения ЗЛП симплекс-методом Рассмотрим задачу из предыдущей лабораторной работы: Z = x1 – 2x2 → min max , 3x1 + 4x2 ≤ 12, 2x1 – x2 ≥ –2, x1 – 2x2 ≤ 6, x1 ≥ 0. При такой постановке эта задача состоит из двух отдельных задач: минимизации Z и максимизации Z. Мы будем рассматривать первую из них – минимизации Z: Z = x1 – 2x2 → min, 3x1 + 4x2 ≤ 12, 2x1 – x2 ≥ –2, x1 – 2x2 ≤ 6, x1 ≥ 0. Сравнив условие задачи и требования при записи канони- ческой формы, видно, что в ней нарушаются все четыре правила. Устраним их. 1) Чтобы от задачи минимизации перейти к задаче на максимум, надо функцию цели умножить на –1. Тогда первая часть исходной задачи в новой постановке будет выглядеть так (–Z) = –x1 + 2x2  max. После решения новой задачи и определения max (–Z) минималь- ное значение Z определяется как min (Z) = –max (–Z). 2) Если ограничение содержит знак ≤, то к его левой части нуж- но добавить некоторую новую неотрицательную переменную. Тогда вместо заданного ограничения мы получаем две записи, включающие собственно ограничение со знаком – и условие неот- рицательности дополнительной переменной. Соответственно, 105 если ограничение имеет знак ≥, то от левой части отнимается новая неотрицательная переменная. Преобразуем ограничения в исходной задаче: 3x1 + 4x2 ≤ 12  3x1 + 4x2 + x3 = 12, где x3 ≥ 0, 2x1 – x2 ≥ –2  2x1 – x2 – x4 = –2, где x4 ≥ 0, x1 – 2x2 ≤ 6  x1 – 2x2 + x5 ≤ 6, где x5 ≥ 0. 3) Если ограничение имеет отрицательный свободный член, то все оно умножается на –1. В данном случае изменит свой вид второе ограничение: –2x1 + x2 + x4 = 2, x4 ≥ 0 4) Если на какую-нибудь переменную не наложено условие неотрицательности, то она заменяется разностью двух неотрица- тельных величин. В исходной задаче условие неотрицательности задано только для переменной x1. Для переменной x2 такое условие не указано. Тогда x2 = (1)2x – (2) 2x , где (1)2x ≥ 0 и (2)2x ≥ 0. Собирая и редактируя снова всю задачу, получаем окончательно (1) (2) 1 2 2( ) 2 2Z x x x     → max, 3x1 + 4 (1)2x – 4 (2) 2x + x3 = 12, –2x1 + (1)2x – (2) 2x + x4 = 2, x1 – 2 (1)2x + 2 (2) 2x + x5 = 6, 1x ≥ 0, (1)2x ≥ 0, (2)2x ≥ 0, x3 ≥ 0, x4 ≥ 0, x5 ≥ 0. Видно, что записав ЗЛП в каноническом виде, изменилась раз- мерность задачи: число неизвестных возросло с двух до шести. 106 8.4.1. Запись начальной жордановой таблицы Чтобы записать начальную жорданову таблицу, преобразуем исходную задачу еще раз, представив функцию цели в виде (–Z) = 1·(–x1) – 2·(– (1)2x ) + 2(– (2) 2x ) + 0·(–x3) +0·(–x4) + 0·(–x5), а ограничения в виде так называемых нуль-строк 0 = 12 + 3·(–x1) + 4·(– (1)2x ) + (–4)·(– (2) 2x ) +1·(–x3), 0 = 2 – 2·(–x1) + 1·(– (1)2x ) –1·(– (2) 2x ) + 1·(–x4), 0 = 6 + 1·(–x1) – 2·(– (1)2x ) + 2·(– (2) 2x ) + 1·(–x5). Нетрудно видеть, что получены они переносом всех слагаемых в правую часть ограничений. Начальная жорданава таблица (табл. 8.1) состоит из верха, где расположены все шесть переменных, трех нуль-строк, соответству- ющих трем ограничениям, и Z-строки. В левом столбце помещены нули и наименование функции цели. Во втором столбце – свобод- ные члены ограничений и функции цели (последний в начальной таблице обычно равен нулю). Условия неотрицатель- ности переменных в таблицу не заносятся. Таблица 8.1 Начальная жорданова таблица –x1 – (1)2x – (2) 2x –x3 –x4 –x5 0= 12 3 4 -4 1 0 0 0= 2 -2 1 -1 0 1 0 0= 6 1 -2 2 0 0 1 (–Z)= 0 1 -2 2 0 0 0 107 8.4.2. Нахождение начального плана Начальный план – это координаты любой вершины ОДР. Чтобы его найти, надо все нули в левом столбце заместить переменными с верха таблицы. Делается это с помощью жордановых исключений. Чтобы некоторая переменная могла быть перенесена, в ее столб- це (исключая элемент Z-строки) должен быть положительный коэффициент; он принимается за разрешающий элемент. Если в столбце два и более положительных коэффициентов, то за разре- шающий элемент берется тот, который соответствует минимально- му отношению свободного члена к коэффициенту столбца. И еще: место в левом столбце, куда собираются перенести неизвестную, не должно быть занято ранее перенесенной неизвестной. При переносе переменной с верха таблицы на ее место должен стать нуль. Это означает умножение столбца на нуль. Поэтому в новой таблице этот столбец вычеркивается. Элементы остальных столбцов и Z-строки преобразуются по второму и четвертому пра- вилам жордановых исключений (см. п. 8.2). В нашей задаче (табл. 8.1) переменная x1 может быть перене- сена в первую или третью строки (надо уточнить, в какую именно), (1) 2x – в первую либо вторую, (2)2x – только в третью, x3 – в первую, x4 – во вторую, x5 – в третью. Начнем с x1. В столбце этой переменной два положительных коэффициента: 3 и 1. Один из них должен быть взят за разрешаю- щий элемент. Но какой именно? Разделим свободный член первой строки 12 на 3, получим 4. А частное от деления свободного члена третьей строки 6 на 1 равно 6, т. е. больше. Поэтому за разрешаю- щий элемент берем 3 (в табл. 8.1 он выделен). Это значит, что x1 переносится в первую строку. Первая жорданова таблица (табл. 8.2) имеет на один столбец меньше. Коэффициенты первой (разрешающей) строки получены по второму правилу жордановых исключений. Коэффициенты вто- рой, третьей и Z-строки преобразованы, как прочие элементы (правило 4). При вычислении коэффициентов первой жордановой таблицы использованы (для наглядности) простые дроби. 108 Таблица 8.2 Первая жорданова таблица – (1)2x – (2) 2x –x3 –x4 –x5 x1 = 4 4/3 -4/3 1/3 0 0 0= 10 11/3 -11/3 2/3 1 0 0= 2 -10/3 10/3 -1/3 0 1 (–Z)= -4 -10/3 10/3 -1/3 0 0 Далее попробуем перенести (2)2x . В соответствующем столбце один положительный коэффициент (10/3), и место в третьей строке не занято. Следовательно, такой перенос возможен. Используя разрешающий элемент 10/3 (в табл. 8.2 он выделен), строим вторую жорданову таблицу (табл. 8.3). Таблица 8.3 Вторая жорданова таблица. – (1)2x –x3 –x4 –x5 x1 = 4,8 0 0,2 0 0,4 0= 12,2 0 0,3 1 1,1 (2) 2x = 0,6 –1 –0,1 0 0,3 (–Z)= –6 0 0 0 –1 Остается заместить нуль во второй строке. Проще всего это сделать, перенеся сюда x4: в соответствующем столбце – один по- ложительный коэффициент и именно во второй строке. Принимаем его за разрешающий и строим третью жорданову таблицу (табл. 8.4). 109 Таблица 8.4 Третья жорданова таблица – (1)2x –x3 –x5 x1 = 4,8 0 0,2 0,4 x4 = 12,2 0 0,3 1,1 (2) 2x = 0,6 –1 –0,1 0,3 (–Z)= –6 0 0 –1 Итак, все нули в левом столбце замещены. Мы пришли в одну из вершин ОДР, называемую начальной. Чтобы найти ее координаты (начальный план), надо перенесенные в левый столбец переменные приравнять соответствующим свободным членам ограничений (x1 = 4,8; x4 = 12,2; (2)2x = 0,6), а оставшиеся на верху таблицы свободные переменные приравнять нулю ( (1)2x = x3 = x5 = 0). Функция цели в этой вершине (–Z = –6). Следует отметить, что если бы мы изменили принятый выше по- рядок переноса неизвестных, то могли бы прийти в другую вершину ОДР. 8.4.3. Нахождение оптимального плана Достигнутую вершину или план x1 = 4,8; (1)2x = 0; (2)2x = 0,6; x3 = 0; x4 = 12,2; x5 = 0, при котором (–Z = –6), надо исследовать на оптимальность. Делается это очень просто. Надо просмотреть коэффициенты Z – строки (в табл. 8.4 это 0; 0; –1). Если среди них имеется хотя бы один отрицательный (а у нас он есть и равен –1), то достигнутая вершина и соответствующий ей план не является оптимальным. Нужно перейти к другой вершине. В качестве разрешающе- го столбца берется тот столбец, в котором находится отрицатель- ный коэффициент функции цели. У нас это последний столбец. 110 В нем три положительных коэффициента. (Если бы в этом столбце не оказалось ни одного положительного коэффициента, то это сви- детельствовало бы о том, что данная задача решения не имеет). Из трех положительных коэффициентов один надо выбрать в каче- стве разрешающего элемента. Для этого находим три частных: 4,8/0,4 = 12; 12,2/1,1 = 11,0909; 0,6/0,3 = 2. Наименьшее из них – последнее. Отсюда разрешающая строка – третья. На ее пересече- нии с разрешающим столбцом стоит разрешающий элемент 0,3, вы- деленный в табл. 8.4. Переходим к четвертой жордановой таблице (табл. 8.5), в кото- рой неизвестная x5 перешла в левый столбец, а на ее место стала (2) 2x . Коэффициенты новой таблицы находятся по четырем прави- лам жордановых исключений. Таблица 8.5 Четвертая жорданова таблица – (1)2x –x3 – (2) 2x x1 = 4 4/3 1/3 –4/3 x4 = 10 11/3 2/3 –11/3 x5 = 2 –10/3 –1/3 10/3 (–Z)= –4 –10/3 –1/3 10/3 Достигли новой вершины с координатами x1 = 4; (1)2x = 0; (2)2x = 0; x3 = 0; x4 = 10; x5 = 2. Функция цели (–Z) выросла с –6 до –4. Этот план также не является оптимальным, так как среди коэффициентов Z – строки имеются отрицательные (–10/3 и –1/3). Надо перейти к следующей вершине ОДР. Поскольку в Z – стро- ке два отрицательных коэффициента, то за разрешающий столбец берем тот, в котором стоит отрицательный коэффициент с большим модулем (–10/3). В разрешающем столбце имеются два положи- тельных коэффициента. Находим частные 4·3/4 = 3 и 10·3/11 = 2,72727. Вторая величина меньше. Поэтому разрешаю- щая строка – вторая, разрешающий элемент равен 11/3. 111 Переходим к пятой таблице (табл. 8.6), меняя местами x4 и (1)2x . Таблица 8.6 Пятая жорданова таблица –x4 –x3 – (2)2x x1 = 4/11 –4/3 1/11 0 (1) 2x = 30/11 3/11 2/11 –1 x5 = 122/11 10/11 3/11 0 (–Z)= 56/11 10/11 3/11 0 Достигли вершины с координатами x1 = 4/11; (1)2x = 30/11; (2) 2x = 0; x3 = 0; x4 = 0; x5 = 122/11. Функция цели выросла до (–Z) =56/11. Т.к. среди коэффициентов Z – строки нет отрица- тельных, полученный план является оптимальным. Переходя к исходным неизвестным, получаем окончательный результат x1 = 4/11 = 0,363636; x2 = (1)2x – (2) 2x = 30/11 = 2,72727, min Z = –max (–Z) = –56/11 = –5,09091. 8.5. Выполнение индивидуальных заданий Используя исходные данные к лабораторной работе № 7, решить свой вариант задания, применяя симплекс-метод. Указания к работе. Функция цели в задаче исследуется на максимум. Условия неот- рицательности на переменные x1 и x2 не накладываются. 112 Литература 1. Бахвалов, Н.С. Численные методы (анализ, алгебра, обыкно- венные дифференциальные уравнения) / Н.С. Бахвалов. – М.: Наука, 1975. – 630 с. 2. Турчак, Л.И. Основы численных методов / Л.И. Турчак. – М.: Наука, 1987. – 320 с. 3. Крылов, В.И., Вычислительные методы: в 2-х томах / В.И. Крылов, В.В. Бобков, П.И. Монастырный. – М.: Наука. – Т.1, 1976. – 333 с. – Т.2, 1977. – 399 с. 4. Караманский, Т. Д. Численные методы строительной механи- ки / Т. Д. Караманский. – М.: Стройиздат, 1981. – 428 с. 5. Амосов, А.А. Вычислительные методы для инженеров: учебное пособие / А. А. Амосов, Ю А Дубинский, Н. В. Копче- нова – М.: Высш. школа, 1994. – 544 с.: ил. 6. Костевич, Л.С. Математическое программирование: Информ. технологии оптимальных решений: учеб. пособие / Л.С. Костевич. Минск: Новое знание. – 2003. – 424 с.: ил. 7. Кузнецов, А.В., Математическое программирование / А.В. Кузнецов, Н.И. Холод. – Минск: Высшая школа, 1984. – 221 с. 8. Мудров, А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль / А.Е. Мудров. – Томск: МП «РАСКО». – 1991. – 272 с.: ил. 9. Дьяконов, В.П. Mathematica 5/6/7. Полное руководство / В.П. Дьяконов. – М.: ДМК, Пресс, 2009. – 624 с. 10. Половко, А.М. Mathematica для студента / А.М. Половко. – СПб.: БХВ-Петербург. – 2007. – 368 с.: ил. 113 ПРИЛОЖЕНИЕ Общие сведения о матрицах Матрицей размера mn называется прямоугольная таблица чисел, расположенных в m строках и n столбцах 11 12 1 21 22 2 1 2 ... ... ... ... n n m n m m mn a a a a a a A a a a           A . Числа, составляющие матрицу, называются элементами матрицы. Матрица размера 1n называется вектором-строкой, а матрица размера m1 – вектором-столбцом. Если m = n (число строк = числу столбцов), то такая матрица называется квадратной матрицей, а n – порядком матрицы. Совокупность элементов квадратной матрицы, расположенных на отрезке, соединяющем левый верхний угол с правым нижним, называется главной диагональю матрицы. Квадратная матрица, у которой все элементы вне главной диаго- нали равны нулю, называется диагональной. Квадратная матрица n-го порядка, у которой на главной диагона- ли стоят единицы, а все остальные элементы равны нулю, 1 0 0 0 1 0 0 0 1         E        . называется единичной. Для ее обозначения используется буква E. Матрица называется нулевой, если все ее элементы равны 0 114 0 0 0 0 0 0 0 0 0         O        . Подобно системам уравнений, квадратные матрицы могут быть симметричными (элементы расположены симметрично относи- тельно главной диагонали), треугольными (элементы, равные нулю, расположены ниже (выше) главной диагонали), ленточными (ненулевые элементы образуют «ленту», параллельную диагонали) и т. д. Определителем n-го порядка называется число ∆, которое символически записывается с помощью квадратной таблицы 11 12 1 21 22 2 1 2 ... ... det ... ... n n m m mn a a a a a a a a a   A и вычисляется следующим образом n = 1, ∆ = |a11| = a11, n = 2, 11 12 11 22 12 21 21 22 a a a a a a a a     , n = 3; 11 12 13 22 23 21 23 21 22 21 22 23 11 12 13 32 33 31 33 31 32 31 32 33 a a a a a a a a a a a a a a a a a a a a a a a a      . Суммой (разностью) двух матриц А и В одинаковых размеров называется третья матрица С тех же размеров, элементы которой определяются по формуле ik ik ikc a b  . 115 Транспонирование матрицы – это замена строк столбцами и наоборот. В публикациях операция транспонирования использу- ется для более удобной записи векторов. Например, вместо того, чтобы изображать вектор в виде 1 2 3 x         удобнее представить его как  1 2 3 Tx  . При умножении (или делении) произвольной матрицы на чис- ло, отличное от нуля, получается матрица тех же размеров, все эле- менты которой получены умножением (или делением) элементов исходной матрицы на это число. Две матрицы А и В можно перемножить между собой и полу- чить третью матрицу С в том и только в том случае, когда число столбцов матрицы А равно числу строк матрицы В. Результирую- щая матрица имеет число строк, как у матрицы А, а число стол- бцов – как у В: m n n p m p   A B C . Элементы результирующей матрицы С определяются по форму- ле 1 n ik ij jk j c a b    , т. е. элементы i-й строки матрицы А перемножа- ются на соответствующие элементы k-го столбца матрицы В и все результаты суммируются. Матричное произведение в общем случае не обладает перестано- вочным свойством. Так, если матрицы А и В из только что рассмот- ренного примера перемножить в обратном порядке В·А, то резуль- тирующая матрица D будет совершенно иной, чем С. Очень часто произведение двух матриц вообще невозможно либо возможно только в одном порядке. Если произведение двух матриц возможно как в прямом, так и в обратном порядке и дает один и тот же ре- зультат, то такие матрицы называются перестановочными. Приме- ром перестановочной матрицы является единичная матрица. Операция перемножения матриц позволяет ввести матричную запись системы линейных уравнений (СЛАУ) 116 11 12 1 1 1 21 22 2 2 2 1 2 ... ... ... ... ... ... n n n n nn n n a a a x b a a a x b a a a x b                                 сокращенно в виде x bA  или AX B , где А – матрица коэффици- ентов; x (X) – вектор (матрица) неизвестных; b (B) – вектор (матрица) свободных членов. Матрицей, обратной матрице A, называется такая матрица А–1, которая удовлетворяет соотношениям 1 1 n n n n n n n n n n          A A A A E . Для того, чтобы обратить квадратную матрицу, она должна быть невырожденной (определитель не равнялся нулю). Невырожденная матрица A является плохо обусловленной, если ее определитель относительно мал, т. е. det A ≈ 0. Использование таких матриц в приближенных вычислениях приводит, как правило, к получению неточных и неустойчивых результатов. 117 Учебное издание ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ СТРОИТЕЛЬСТВА Лабораторный практикум для студентов специальности 1-70 02 01 «Промышленное и гражданское строительство» Составители: СТРЕЛЮХИН Александр Владиславович БОГОМОЛОВА Гелена Станиславовна СОРОКИНА Елена Леонтьевна Редактор О. В. Ткачук Компьютерная верстка Ю. С. Кругловой, А. Е. Дарвина Подписано в печать 31.05.2016. Формат 6084 1/16. Бумага офсетная. Ризография. Усл. печ. л. 6,74. Уч.-изд. л. 5,27. Тираж 100. Заказ 468. Издатель и полиграфическое исполнение: Белорусский национальный технический университет. Свидетельство о государственной регистрации издателя, изготовителя, распространителя печатных изданий № 1/173 от 12.02.2014. Пр. Независимости, 65. 220013, г. Минск.