Прикладная математика.

© Лекции подготовлены доц. Мусиной М.В. Прикладная математика. Наступивший век - это век всеобщей информатизации. Математические расчеты, проводимые с помощью компьютеров, прочно проникли в самые разнообразные научные дисциплины, технику, экономику, управление и другие сферы человеческой деятельности, причем не только в традиционные области приложения математики, но и туда, где математика используется совсем недавно. Применяемый математический аппарат стал намного разнообразнее. Если еще 30 – 50 лет назад инженерам требовались знания, лежащие в пределах школьной математики, дифференциального и интегрального исчислений, то теперь применяется аппарат ранее неизвестный и даже целые математические дисциплины, например линейное программирование, созданные в первую очередь с целью приложения. В результате повысились требования к математической подготовке инженера. Это связано с все большим распространением компьютерного (модельного) подхода к решению инженерных задач. Понятие «моделирование» стало за последнее десятилетие едва ли не самым распространенным в естественно - научной и технической литературе. Сегодня трудно представить проектную или конструкторскую организацию, не использующую моделирования в своей практике в той или иной мере. И чем сложнее проектируемый объект, тем важнее роль моделирования в его изучении и создании. Нет никакого сомнения в том, что развитие и применение математических моделей и математического аппарата будут только еще более усиливаться. Этим объясняется возросший интерес к вопросу о том, как создаются математические модели, как они изучаются и как интерпретируются. Моделирование – один из основных методов познания, является формой отражения действительности и заключается в выяснении или воспроизведении тех или иных свойств реальных объектов предметов и явлений с помощью других объектов, процессов, явлений, либо в виде изображения, совокупности уравнений, алгоритмов и программ. Моделью называется объект-заместитель, который в определенных условиях может заменять объект-оригинал, воспроизводя интересующие нас свойства и характеристики оригинала, причем имеет существенные преимущества: - дешевизну; - наглядность; - легкость оперирования и т.п. Модель – образ реальной системы (объекта, процесса) в материальной или теоретической форме. Этот образ отражает существенные свойства объекта, он замещает реальный объект в ходе исследования и управления. Модель отображает (описывает, воспроизводит) некоторые интересующие черты исследуемого объекта. Типы моделей: физические или аналоговые (модели самолетов в аэродинамической трубе, план, карта, график) и математические модели, где количественно определенные переменные связаны математическими законами ( модели в физике, экономике). Своего рода графической моделью является проект землеустройства, топологическая карта, схема и пр. Математическая модель- это приближенное описание какого-либо класса явлений внешнего мира, выраженное с помощью математической символики, то есть система математических соотношений, приближенно в абстрактной форме описывающих изучаемый процесс или систему. Создание математической модели. При построении модели реальное явление неизбежно упрощается, схематизируется, и эта схема описывается с помощью того или другого математического аппарата. Чем удачнее будет подобрана математическая модель, чем лучше она будет отражать характер- © Лекции подготовлены доц. Мусиной М.В. ные черты явления, тем успешнее будет исследование и полезнее – вытекающие из него рекомендации. Общих способов построения математических моделей не существует. В каждом конкретном случае модель выбирается исходя из задачи исследования (какие параметры требуется определить и влияние каких факторов отразить). Математическая модель должна отражать важнейшие черты явления и вместе с тем должна быть по возможности простой. Основные этапы создания математической модели. При создании модели встает вопрос о решении математической задачи. В какой бы математической форме не была реализована модель в дифференциальной, линейной или другой, необходимо найти решение. В качестве компонентов математических моделей часто встречаются такие как решение различного вида уравнений и их систем, а также вычисление интегралов. И далеко не все эти математические задачи можно решить аналитически. Здесь и вступает в свои права прикладная математика. Раздел математики, имеющий дело с созданием и обоснованием численных алгоритмов для решения сложных задач различных областей науки, называют прикладной математикой. Главная задача прикладной (вычислительной) математики – фактическое нахождение решения с требуемой (оцениваемой) точностью. Этим она отличается от классической математики, которая основное внимание уделяет исследованию условий существования и свойств решения. В истории прикладной математики можно условно выделить 3 периода. Первый начался 3 – 5 тысячи лет назад. Он был связан с ведением конторских книг, вычислением площадей и объемов, расчетами простейших механизмов; иными словами – с несложными задачами арифметики, алгебры и геометрии. Вычислительными средствами служили сначала собственные пальцы, а затем – счеты. Исходные данные содержали мало цифр, и большинство выкладок выполнялось точно, без округлений. Второй период начался с Ньютона. В этот период решались задачи астрономии, геодезии и расчета механических конструкций, сводящиеся либо к обыкновенным дифференциальным уравнениям, либо к алгебраическим системам с большим числом неизвестных. Вычисления выполнялись с округлением; нередко от результата требовалась высокая точность, так что приходилось сохранять до 8 значащих цифр. © Лекции подготовлены доц. Мусиной М.В. Вычислительные средства стали разнообразнее: таблицы элементарных функций, затем арифмометр и логарифмическая линейка. Скорость этих средств была невелика, и вычисления занимали дни, недели и даже месяцы. Третий период начался примерно в 40-х годах 20 века. Военные задачи требовали недоступных человеку скоростей и привели к разработке ЭВМ. Скорость даже простейших ЭВМ настолько превосходила скорость механических средств, что стало возможным проводить вычисления огромного объема. Это позволило численно решать новые классы задач. Для решения многих задач до сих пор используются методы разработанные в «доэлектронный» период, но также появилось множество новых численных методов. Целью изучения данной дисциплины являются:  Знакомство с принципами построения алгоритмов и методикой постановки задач для приближенного решения на ЭВМ.  Изучение вычислительных методов, применяемых при решении прикладных задач, не имеющих аналитического решения, либо имеющих его, но, по ряду причин, получение которого затруднено;  Ознакомление с основными источниками погрешностей, их оценкой и методами устранения; Сложные вычислительные задачи, которые необходимы для получения решения математических моделей, можно разбить на ряд элементарных – таких как вычисление интеграла, решение дифференциального уравнения, решение нелинейных уравнений и систем. Часть из них может быть решена аналитически, однако для большинства используются приближенные вычислительные методы. В курсе высшей математики мы уже познакомились с таким приближенным методом как разложение функций в степенные ряды Тейлора. Распространенным способом решения является применение численных алгоритмов (численных методов). Их существует великое множество и первоначально они не были связаны с применением ЭВМ, а осуществлялись также вручную. Появление компьютеров сделало очень удобным и простым (поскольку они используются и в пакетах математических программ) их применение. В настоящее время самые распространенные методы записаны в книгах уже в виде готовых программ, «численных рецептов». Мы воспользуемся этими рецептами для некоторых распространенных задач: решения нелинейных уравнений, систем линейных уравнений, вычисления определенного интеграла, решения дифференциальных уравнений, методам численной аппроксимации функций. Макросы VBA Excel и их базовые элементы. Следует сказать, что разработанный численный метод содержит обычно только принципиальную схему решения задачи, не включая детали, без которых невозможна реализация метода на ЭВМ. Необходимо подробно детализировать все этапы вычислений, то есть составить программу на выбранном языке программирования, например Фортран, Си и пр. Мы будем реализовать решение вычислительных задач в виде макросов на языке VBA (Visual Basic for Application). Программа представляется в виде последовательности строк, описывающих процесс (алгоритм) решения поставленной задачи. Каждая из программных строк состоит из одного или нескольких операторов. Макросом (или макрокомандой) в EXCEL называется последовательность операций, снабженная именем, которая может быть выполнена автоматически. Программа, написанная на VBA EXCEL – макрос. Кроме этого можно называть такую программу приложением или проектом. Для того чтобы пользоваться программами (макросами) не обязательно быть программистом, тем более что при изучении дисциплины будут использоваться готовые программы. Однако, даже при использовании готовых макросов, надо понимать, что записано в тексте кода(текста программы), какая последовательность действий будет осуществляться, в каких строках и почему необходимо сделать изменения. Поэтому мы познакомимся с базовыми конструкциями, имеющимися во всех алгоритмических языках – операторами © Лекции подготовлены доц. Мусиной М.В. присваивания, условными операторами, циклами, массивами и подпрограммами, с тем, чтобы была возможность модифицировать макросы для конкретных нужд. Операторы ввода – вывода, присваивания. Пример простейшей программы: Sub Пифагор() a=3 b=4 c = Sqr(a ^ 2 + b ^ 2) End Sub Первая строка кода заглавие программы, с этой строки начинается любая программа. Здесь содержится имя программы, в данном случае Пифагор(). Любая программа заканчивается оператором End Sub. В блок – схеме алгоритма этим двум строкам соответствует обозначение: Внутри этих двух строк операторы программы. Мы записали 3 оператора которые называются операторы присваивания. Общий вид: a = b . Справа от равенства записывается алгебраическое выражение, вычисленное значение которого будет присвоено переменной стоящей слева. Арифметические действия в выражениях совпадают с аналогичными операциями в EXCEL. Математическая запись одного из операторов c = √a2 + b2. a, b, c – переменные. Переменные в программе имеют такой же смысл, как и в математике. Перед тем как использовать переменную ее рекомендуется описать (объявить). Синтаксис оператора описания переменной: Dim переменная [As тип]. Здесь Dim – ключевое слово описывающее переменную. (сокращение от слова dimension – размер). Тип данных может быть целым (Integer), вещественным (Single), текстом (String) и пр. В простейших программах строго описывать необходимо только массивы. Кроме математических действий в операторе присваивания могут содержаться математические функции. Например: c = Sqr(a ^ 2 + b ^ 2) эквивалентно математической записи c  a2  b 2 . Некоторые математические функции, используемые в программах: sin(x), cos(x), sqr(x) ( x ), exp(x) (ex), log(x) (ln(x)), abs(x) ( x ), atn(x) (arctg(x)). Часть математических функций в VBA отсутствует и их необходимо рассчитывать по формулам: arcsin x   atn x / sqr 1  x * x   arccos x    atn x / sqr 1  x * x  2 chx   exp x   exp  x  / 2 shx   exp x   exp  x  / 2 Исполняемые операторы в блок – схемах обозначаются: В примере простейшей программы двум переменным a и b присваиваются конкретные значения. Чтобы рассчитать переменную с уже с другими значениями нужно внести изменения в программу. Это неудобно. Действия в программе должны выполняться многократно с разными значениями входных данных. Поэтому значения переменных a и b вводятся в программу с листа MS EXCEL. © Лекции подготовлены доц. Мусиной М.В. a = Range("a").Value b = Range("b").Value Вывод результатов осуществляется на лист MS EXCEL. Range("c6").Value = c В блок – схемах алгоритмов операции ввода – вывода соответсвует диаграмма: Условные операторы и операторы перехода. Операторы рассмотренной программы выполнялись поочередно. Это линейная программа. Для изменения последовательности выполнения операторов (т.е. для ветвления) программы используются операторы перехода. Оператор безусловного перехода: GOTO метка («номер строки») Если в определенном месте программы необходимо выполнить те или иные операторы в зависимости от некоторых условий, то применяются операторы условного перехода. IF условие THEN оператор Под условием понимается некоторое логическое выражение (высказывание), которое может быть либо истинным (TRUE – логическая 1), либо ложным (FALSE – логический 0). Условие часто записывается в виде арифметического отношения с использованием конструкций “>”,”<”,”>=”,”<=”, “< >”(). Если условие истинно, то оператор выполняется, если ложно, то выполняется следующий. Пример: IF (X > 12 Or X < 5) THEN X = X + 1 Если при истинности условия требуется выполнить несколько операторов, то используют конструкцию вида: IF условие THEN операторы END IF В случае, когда имеется альтернатива, то условный оператор приобретает вид: IF условие THEN операторы 1 ELSE операторы 2 END IF В блок-схемах условные операторы обозначаются: Оператор цикла и пользовательские функции. Для многократного выполнения блока операторов необходимо использовать оператор цикла. Когда количество выполнений заданного блока операторов известно заранее то используется оператор FOR…TO…STEP…NEXT. Вид оператора: FOR счетчик = начало TO конец [STEP шаг] операторы NEXT Счетчик это переменная, которой присваивается значение начало. Затем значение этой переменной увеличивается на величину шага (по умолчанию шаг равен 1). Операторы цикла выполняются до тех пор, пока значение счетчика не достигнет конечного значения. © Лекции подготовлены доц. Мусиной М.В. Пример: For i = 1 To 10 x=a+i·h Next i Пользовательские функции – это блок операторов, предназначенный для многократного выполнения в разных точках программы, оформляется как новая процедура (программа). Процедуры делятся на функции и подпрограммы. Пользовательская функция: FUNCTION название (параметры) [As тип] операторы END FUNCTION Пользовательская подпрограмма: SUB название (параметры) операторы END SUB Обращение к подпрограмме: CALL название (параметры) Пример: Function fnf(x) As Single fnf = Sqr(2 * x + 1) End Function Решение нелинейных уравнений. Постановка задачи: Найти коэффициент погрешности прибора σ при проведении геодезических измерений из уравнения: δ · cos σ – υ · σ 2 + η = 0 Значения δ = 0,186, υ = 4,18, η = 0,124. Т.е необходимо решить уравнение: 0,186 · cos σ – 4,18 · σ 2 + 0,124 = 0 Это нелинейное уравнение, для которого нет аналитического метода. Для его решения можно применить численные методы математики. Общие понятия. Всякое уравнение с одним неизвестным можно записать в виде f(x) = 0 (1) Корнем (или решением) уравнения (1) называется значение x такое, что f(x) = 0. Уравнения делятся на алгебраические и трансцендентные (т.е. уравнения, содержащие тригонометрические функции или другие специальные функции типа ex, ln x и т.д.). Пример алгебраич. уравнения: xn + an-1xn-1 + … + a1x + a0 = 0. (2) Для уравнения (2) явные формулы найдены только до n = 4. Для уравнений более высоких степеней точных решений не существует (т. Абеля). Т.е. под задачей отыскания решений уравнения (1) будем понимать задачу вычисления с заданной точностью e конечного числа подлежащих определению корней этого уравнения (обычно задачу уточняют, т.е. какие именно корни нужно найти – положительные, максимальные, или корни из заданного интервала). Корни могут быть простые и кратные. Методы решения прямые (с помощью формул, позволяющих найти точное решение) и итерационные (т.е. многократное повторение некоторого алгоритма с нахождение определенного приближенного решения. Основные этапы решения: 1) Локализация (отделение) корней. 2) Уточнение решения (итерационное уточнение корня) Первый этап – локализация корней. © Лекции подготовлены доц. Мусиной М.В. Определение. Отрезок [a, b] содержащий только один корень x уравнения (1) называется отрезком локализации корня x. Т.е. под отделением корня понимают нахождение отрезка, на котором лежит этот и только этот корень данного уравнения. Способы отделения корней многообразны. Иногда из физ. сообр. графический метод. Построение таблиц функции. При этом о локализации корня судим по перемене знаков. Несколько теорем из математического анализа. Теорема 1. Если функция y = f(x) непрерывна на отрезке [a, b] и принимает на его концах значение разных знаков, т.е. f(a) f(b) < 0, то внутри отрезка [a, b] существует, по крайней мере, один корень уравнения (1). Теорема 2. Если функция y = f(x) непрерывна на отрезке [a, b], f(a) f(b) < 0 и f '(x) на интервале сохраняет знак, то внутри отрезка [a, b] существует единственный корень уравнения f(x) = 0. Примеры локализации корней. Графический: 4(1 – x2) – ex = 0 => 1 – x2 = 0,25 ex Строим график: (рисунок) Отрезки [-1, 0] и [0, 1] являются отрезками локализации. Пример. Локализация корней уравнения x3 – 1,1x2 – 2,2x + 1,8 = 0 с помощью таблицы на отрезке [-2, 2] x f(x) -2 -6,2 -1,6 -1,592 -1,2 1,128 -0,8 2,344 -0,4 2,44 1,8 0,4 0,808 0,8 -0,152 1,2 -0,696 Функция меняет знак на концах отрезков [-1,6; -1,2]; [0,4; 0,8]; [-1,6; 2,0] Второй этап. Итерационное уточнение корня. Если найдены отрезок [a, b] на котором есть один корень уравнения, то любую точку отрезка можно принимать за приближенное значение корня. При этом погрешность приближения не превышает длины [a, b]. Т.е. задача нахождения приближенного значения сводится к нахождению отрезка [a, b] ( b  a   ) содержащий только один корень уравнения. 1. Метод половинного деления (бисекции). За приближенное значение можно принять середину отрезка [a, b]. c=(a+b)/2 . Проверяя знаки f(a), f(b), f(c) выясним в каком из отрезков [a, c] или [c, b] содержится корень x*  [a, c] , если f(a)f(c) < 0 ; x*  [c, b] , если f(c)f(b)<0 . Выбранный отрезок принимаем за [a, b] и повторяем это до тех пор, пока получаемый отрезок не сожмется до заданной степени точности ( b  a   ). За приближенное значение можно принять середину отрезка [a, b] . При n итерациях получим соотношение (b - a)/2n   из которого можно вычислить число итераций, необходимое для достижения заданной степени точности n ln2(b - a)/ . Геометрическая интерпретация. Решим нелинейное уравнение: 0,186 · cos σ – 4,18 · σ 2 + 0,124 = 0 , заменив переменную σ на x. Т.е. f(x) = ,186 · cos x – 4,18 · x 2 + 0,124 Для локализации корня строим таблицу и график. Очевидно, что корень находится на отрезке [0,2; 0,3] 1 -0, © Лекции подготовлены доц. Мусиной М.В. 0,4 0,2 -0,3 -0,2 -0,1 0,1 -0,2 -0,4 -0,6 -0,8 -1 Программа (макрос): Sub bis() Dim a As Single, b As Single, eps As Single a = Range("a").Value b = Range("b").Value eps = Range("eps").Value 10 fa = fnf(a) fb = fnf(b) c = (a + b) / 2! fc = fnf(c) If fc = 0 Then GoTo 20 If Abs(a - b) < eps Then GoTo 20 If fa * fc < 0 Then b = c Else a = c GoTo 10 20 Range("c6").Value = c Range("fc").Value = fc End Sub Function fnf(x) As Single fnf = 0.186 * Cos(x) - 4.18 * x * x + 0.124 End Function 0,2 0,3 0,4 0,5 0,6 © Лекции подготовлены доц. Мусиной М.В. Результат расчета: x = 0.269 Т.е. значение коэффициента σ = 0,269. 2. Метод простой итерации (последовательных приближений). Этот метод является наиболее общим и многие другие методы можно представить как некоторую вариацию метода простых итераций. Для применения этого метода уравнение f(x) = 0 (1) преобразуется к виду: x = φ(x). (3) Это можно сделать, различными способами, например, x =x +c· f(x), с  0. φ(x) – итерационная функция. Предположим, что выбрано некоторое начальное приближение корня уравнения (1) x0 (x0  [a,b]). Подставим его в правую часть уравнения (3). Тогда получим некоторое число x1 = φ(x0) Подставляя теперь в правую часть равенства (1) вместо x0 число х1, получим новое число x2 = φ(x1) Повторяя этот процесс, будем иметь последовательность чисел . xn = φ(xn-1) , n = 0, 1, 2… (4) Такая последовательность называется итерационной. Если эта последовательность - сходящаяся, т. е. существует предел x *  lim x n , то n  этот предел является корнем уравнения и может быть вычислен по формуле (3) с любой степенью точности. Графическая интерпретация метода: решением уравнения (3) будет абсцисса точки пересечения прямой y=x с кривой y= φ(x). © Лекции подготовлены доц. Мусиной М.В. При применении данного метода встает вопрос об условиях его сходимости. Теорема. Пусть функция φ(х) определена и непрерывно дифференцируема на некотором отрезке [а,b] и ее производная |φ΄(x)| ≤ q, где 0 ≤ q<1 , тогда итерационная последовательность сходится к единственному на [а,b] корню уравнения. Скорость сходимости метода зависит от величины q. Чем она меньше, тем быстрее сходится метод. Поэтому на практике при нахождении корней методом простой итерации стараются таким образом подобрать параметр c в формуле, чтобы производная φ΄(x) в окрестности корня по абсолютной величине была как можно меньше. Критерием окончания расчета при заданной точности ε можно считать выполнение неравенства | xn - xn-1 | < ε. Пример 1. Найти корни уравнения ex – 10x =0 с точностью ε = 0,0001. Корни отделяем графически. Корни находятся на отрезках [0, 1] и [3, 4]. 60 50 40 30 20 10 -1 -0,5 0,5 1 1,5 2 2,5 3 3,5 4 4,5 -10 Преобразуем уравнение к виду: x = 1/10·ex Условия сходимости выполняются на отрезке [0, 1]. За начальное приближение примем x0 = 1. © Лекции подготовлены доц. Мусиной М.В. Запишем результаты итераций в таблицу. xn | xn - xn-1 | N итерации 1 2 3 4 5 6 1 0,271828 0,131236 0,114024 0,112078 0,11186 0,111836 0,72817 0,14059 0,01721 0,00195 0,00022 2,4E-05 Корень уравнения x = 0,11183. Пример 2. Решим уравнение ex·sinx – 1 =0. Отрезок локализации корня 0;  2  . Для того чтобы подобрать итерационную функцию будем представлять уравнение в виде x =x +c· f(x). x  x  c e x  sin x  1  x   x  c e x  sin x  1 Для того чтобы условие сходимости выполнялось, необходимо таким образом подобрать параметр c, чтобы |φ΄(x)| ≤ 1.  x   1  ce x sin x  cos x   1     Очевидно, что условие не может выполняться, если c > 0. Чтобы понять как ведет себя производная функции φ(х), найдем  x  .  x   2ce x cos x На интервале 0;  2  ex > 0, cosx > 0, c < 0. Т.е  x   0 . Следовательно φ΄(x) убывает. Тогда наибольшее значение функции φ΄(x) на интервале при x = 0: φ΄(0) = 1 + c .    А наименьшее значение при x  :     1  4,8 c . 2 2 1  c  0  c  1 , 1  4,8 c  0  c  0,21 . В качестве c выберем среднее арифметическое двух значений: с = - 0,6.  x   x  0,6 e x  sin x  1 За начальное приближение выберем середину отрезка локализации x0 = 0,7. xn | xn - xn-1 | N итерации  1 2 3 4 5 6 7 8 9 10 11 12 0,7 0,521623 0,617926 0,573088 0,596025 0,584744 0,59041 0,587593 0,589001 0,588299 0,588649 0,588475 0,588562  0,178377 0,096303 0,044838 0,022938 0,011282 0,005667 0,002818 0,001408 0,000702 0,00035 0,000175 8,73E-05 Корень уравнения x = 0,5886. 3. Метод Ньютона (касательных). Зададим некоторое начальное приближение x0 [a,b] и линеаризуем функцию f(x) в окрестности x0 с помощью отрезка ряда Тейлора © Лекции подготовлены доц. Мусиной М.В. f(x) = f(x0) + f '(x0) (x - x0). (5) Вместо уравнения (1) решим линеаризованное уравнение f(x0) + f '(x0) (x - x0) = 0, трактуя его решение x как первое приближение к корню x1 = x0 - f(x0)/f '(x0) . Продолжая этот процесс, приходим к формуле Ньютона xn+1 = xn - f(xn)/f '(xn) (6), которую можно считать итерационным процессом с итерирующей функцией φ(х)= x - f(x)/f '(x). Геометрическая интерпретация этого процесса показана на рис. Уравнение (5) является уравнением линии касательной к кривой f(x) в точке x0, поэтому этот метод называют методом касательных. Формулу (6) можно легко вывести также из геометричеcких соображений. Рис. Если нулевое приближение выбрано достаточно близко к корню, ньютоновские итерации сходятся. Скорость сходимости квадратичная, то есть на каждом шаге погрешность пропорциональна квадрату предыдущей. Следует отметить, что улучшение сходимости этого метода по сравнению с предыдущими достигается увеличением затрат на выполнение каждого шага, так как в каждом шаге требуется вычислять не только значение функции f(x), но и ее производной f '(x). Практический критерий окончания расчета при заданной точности ε - выполнение неравенства | xn - xn-1 | < ε. Для нахождения корней произвольной дифференцируемой функции чаще всего применяют метод Ньютона, особенно если известны разумные начальные приближения для корней. Пример. Рассмотрим решение уравнения f(x) ≡ x2 – a = 0. Формула метода Ньютона принимает вид: x2  a 1  a  xn 1  xn  n   xn   2 xn 2 xn  Мы получили формулу, которая позволяет очень быстро находить квадратный корень только с помощью умножения и деления. Вот, например 5 , найденный по этой формуле. N итерации 1 2 3 4 5 6 7 xn 1 3 2,333333 2,238095 2,236069 2,236068 2,236068 2,236068 © Лекции подготовлены доц. Мусиной М.В. Решение систем линейных уравнений. Пусть дана система n линейных алгебраических уравнений с n неизестными.  a11 x1  a12 x2  ...  a1n xn  b1 a x  a x  ...  a x  b  21 1 22 2 2n n 2 (1)  ...   an x1  an2 x2  ...  ann xn  bn и пусть detA ≠ 0 (определитель матрицы системы ≠ 0). В этом случае система имеет единственное решение. Тогда в матричной форме АX = B. На практике кроме существования и единственности решения важна еще устойчивость относительно погрешностей правой части и элементов матрицы. Иногда небольшое изменение исходных данных может сильно изменить решение. В этом случае систему называют плохо обусловленной. Очевидно, у плохо обусловленных систем det A ≈ 0. Геометрически плохо обусловленная система из 2-х уравнений соответствует почти параллельным прямым. Пример. 5 x  7 y  12 Рассмотрим систему уравнений:  . Эта система имеет единственное 7 x  10 y  17 решение x = 1 и y = 1. Теперь возьмем значения неизвестных x = 2,415 и y = 0. Подставим эти значения в исходную систему. 5 x  7 y  12,075  7 x  10 y  16,905 Если округлить правые части этой системы до двух знаков, то она будет совпадать с исходной. И если исходные данные будут заданы с такой точностью, то это решение так же хорошо отвечает условиям поставленной задачи как и x = 1 и y = 1. На графике видно, что прямые линии заданные уравнениями системы практически параллельны. 2 1,5 1 0,5 -0,5 0,5 1 1,5 2 2,5 3 © Лекции подготовлены доц. Мусиной М.В. Это пример плохо обусловленной системы. В этом случае найти численное решение системы трудно. Методы решения систем делятся на прямые (точные) и итерационные. Прямыми методами называются методы, позволяющие получить решение системы (1) за конечное число арифметических операций. К этим методам относятся метод Крамера, метод Гаусса, LU-метод и т.д. Следует заметить, что реализация прямых методов на компьютере приводит к решению с погрешностью, т.к. все арифметические операции над переменными с плавающей точкой выполняются с округлением. Точные методы применяются если размер системы не очень велик до 103 уравнений. Метод Гаусса. Метод Гаусса – метод последовательных исключений неизвестных, это один из самых распространенных методов решения линейных алгебраических уравнений. Он известен в различных вариантах уже более 2000 лет. Вычисления по методу Гаусса состоят из двух основных этапов: прямого хода и обратного хода. Прямой ход – последовательное исключение неизвестных из системы для преобразования ее к эквивалентной системе с верхней треугольной матрицей. Вычисления значений неизвестных производят на этапе обратного хода. Самый простой вариант метода Гаусса это схема единственного деления. Прямой ход. 1-й шаг. Целью этого шага является исключение неизвестного x1 из уравнений i = 2, 3 … n. Предположим, что а11 ≠ 0. а11 – главный (или ведущий) элемент первой строки. Найдем множители первого шага a mi1   i1 . a11 Последовательно умножая первое уравнение на mi1 и складывая с i-м уравнением, исключим x1 из всех уравнений кроме первого. Получим систему a11 x1  a12 x2  ...  a1n xn  b1 (1 ) a22 x2  ...  a2(1n) xn  b2(1 ) эквивалентную исходной в которой первое неизвест..................................... (1 ) (1) an(12) x2  ...  ann xn  bnn ное исключено из всех уравнений кроме первого. Формулы для расчета новых коэффициентов системы: aij1   aij  mi1a1 j , bi1   bi  mi1b1 (i, j = 2, 3, … n) 2 шаг. Исключаем неизвестное x2 из уравнений с номерами i = 3, 4, …, n. Аналогичным 1  0 . Множители второго шага образом, считаем главным элементом a22 1  mi 2   ai12 a22 (i = 3, 4, … , n) Умножаем второе уравнение на mi2 и складывая с i-м уравнением, исключим x2 из всех уравнений кроме первого. a11 x1  a12 x2  a13 x3  ...  a1n xn  b1 (1 ) (1) a22 x2  a23 x3  ...  a2(1n) xn  b2(1 ) (2 ) a33 x3  ...  a3( 2n) xn  b32  ..................................... (2 ) (2 ) an(23) x3  ...  ann xn  bnn Коэффициенты рассчитываются по формулам: aij2   aij1   mi 2 a2 j , bi2   bi1   mi 2 b22  (i, j = 3, 4, … n) Все остальные шаги выполняются аналогично. Таким образом на k – том шаге расчета новые коэффициенты системы рассчитываются по формулам © Лекции подготовлены доц. Мусиной М.В. k 1  aijk   aijk 1   mik akj , bik   bik 1   mik bkk11  , где mik   aikk 1  akk (i, j = k+1, … n) После n-1 шага вычисления прямого хода заканчиваются и, последовательно исключив все неизвестные, получим систему треугольного вида: a11 x1  a12 x2  ...  a1k xk  ...  a1n xn  b1, (1 ) a22 x2  ...  a2(1k) xk  ...  a2(1n) x n  b2(1 ) ..................................... an( n11,n)1 xn 1  an( n11,n) xn  bn( n11) ( n 1 ) ann xn  bn( n 1) Заметим, что выполнение прямого хода было возможно при условии, что все aii(l ) , l  1,2...m  1 не равны нулю. Если условие не выполняется, то необходимо переставить строки так, чтобы оно выполнялось. Обратный ход. bn( n 1 ) Из последнего уравнения системы легко найти xn. xn  ( n 1) ann Выполняя последовательные подстановки в системе можно получить все значения неизвестных xn-1, xn-2, …, x1 по формулам: n 1 xi  ( i 1 ) (bi( i 1 )   aij( i 1 ) x j ) (i = n-1, n-2, …, 1) aii j  i 1 Эта процедура получила название обратный ход метода Гаусса. Пример.  x1  0,17 x2  0,25 x3  0,54 x4  0,3   0,47x1  x2  0,67 x3  0,32 x4  0,5  - 0,11x1  0,35 x2  x3  0,74 x4  0,7 - 0,55x1  0,43 x2  0,36 x3  x4  0,9  Решение: Запишем в таблицу расширенную матрицу системы, с которой будем производить действия. 1 0,47 -0,11 -0,55 0,17 1 0,35 0,43 0,25 0,67 1 0,36 0,54 -0,32 -0,74 1 0,3 0,5 0,7 0,9 Последовательно умножим первую строку на 0,47, -0,11 и – 0, 55 и отнимем соответственно из 2, 3 и 4 строки. Преобразованная матрица: 1 0,17 0,9201 0,3687 0,5235 0,25 0,5525 1,0275 0,4975 0,54 -0,5738 -0,6806 1,297 0,3 0,359 0,733 1,065 Вторую строку разделим на 0,92 и получим множители второй строки 1 0,600478 -0,62363 0,390175 Далее последовательно умножаем эту строку на 0, 36 и 0,52 и вычитаем из 3-й и 4-й строк. 1 0,17 1 0,25 0,600478 0,806104 0,18315 0,54 -0,62363 -0,45067 1,623469 0,3 0,390175 0,589142 0,860743 Множители третьей строки: 1 -0,55907 0,730852 Эту строку умножаем на 0,183 и вычитаем из 4-й строки. 1 0,17 0,25 0,54 0,3 © Лекции подготовлены доц. Мусиной М.В. 1 0,600478 1 -0,62363 -0,55907 1,725863 0,390175 0,730852 0,726888 Теперь разделим элементы 4-й строки на 1,72 и получаем матрицу треугольного вида: 1 0,17 1 0,25 0,600478 1 0,54 -0,62363 -0,55907 1 0,3 0,390175 0,730852 0,421174 Осуществляем обратный ход метода Гаусса. x1 = 0,181; x2 = 0,0725; x3 = 0,966; x4 = 0,421 Метод Гаусса может быть легко реализован на компьютере. При выполнении вычислений, как правило, не интересуют промежуточные значения матрицы А. Поэтому численная реализация метода сводится к преобразованию элементов массива размерности (n×(n+1)), где n+1 столбец содержит элементы правой части системы. Один из основных недостатков метода Гаусса связан с тем, что при его реализации накапливается вычислительная погрешность. Для того, чтобы уменьшить рост вычислительной погрешности применяются различные модификации метода Гаусса. Например, метод Гаусса с выбором главного элемента по столбцам, в этом случае на каждом этапе прямого хода строки матрицы переставляются таким образом, чтобы диагональный угловой элемент был максимальным. При исключении соответствующего неизвестного из других строк деление будет производиться на наибольший из возможных коэффициентов и следовательно относительная погрешность будет наименьшей. Существует метод Гаусса с выбором главного элемента по всей матрице. В этом случае переставляются не только строки, но и столбцы. Использование модификаций метода Гаусса приводит к усложнению алгоритма увеличению числа операций и соответственно к росту времени счета. Модифицированные методы Гаусса предлагается изучить самостоятельно, как дополнительный материал. Выполняемые в методе Гаусса преобразования прямого хода, приведшие матрицу А системы к треугольному виду позволяют вычислить определитель матрицы a11 a12  a1m (1) 0 a22  a2(1m) (1 ) det A   a11  a22  am( m,m1 ) .     0  am( m,m1) Метод Гаусса также позволяет найти обратную матрицу. Решение систем линейных уравнений методом простых итераций. Если система имеет большую размерность (103 – 106 уравнений) или матрица системы разрежена, более эффективны для решения непрямые итерационные методы. Разреженная матрица возникают в системе, где многие коэффициенты при неизвестных равны нулю. Итерационные методы (методы последовательных приближений) состоят в том, что решение системы (1) находится как предел последовательных приближений x (n ) при n   , где n номер итерации. При использовании методов итерации обычно задается некоторое малое число   0 и вычисления проводятся до тех пор, пока не будет выполнена оценка x ( n)  x   . К этим методам относятся метод Зейделя, Якоби, метод верхних релаксаций и т.д. Рассмотрим один из простых методов – метод простой итерации. © Лекции подготовлены доц. Мусиной М.В. Пусть дана система  a11 x1  a12 x2  ...  a1n xn  b1 a x  a x  ...  a x  b  22 2 2n n 2 вид:  21 1 (1). ...   an x1  an2 x2  ...  ann xn  bn линейных уравнений, имеющая Эта система эквивалентна векторной (матричной) записи A  X  B , где А – матрица системы,  x1     x2  X    - вектор – столбец неизвестных,    x   n  b1     b2  B    - вектор – столбец свободных членов.    b   n Данную систему можно преобразовать к виду X  C  X  F (2), где С – некоторая матрица, F - вектор – столбец. Обозначим C  X  F   X , то есть система принимает форму X   X . То есть, начиная с некоторого начального значения (нулевого приближения) получим последовательность. 1  0  X   X    2  1  X   X    … k 1  k  X   X  .   Аналогичную последовательность мы получали при решении нелинейного уравнения методом простой итерации x = φ(x). Возникает только вопрос о том, как привести систему (1) к виду удобному для итераций (2). В общем случае это не простая задача, требующая специальных знаний, но в некоторых случаях можно поступить очень просто. Например, можно сделать такие преобразования: A  X  A  I  I   X  B , где I – единичная матрица.     A  I   X  I  X  B X  I  A   X  B Самый простой способ следующий. Из первого уравнения системы (1) выразим неизвестное x1: 1 x1  b1  a12 x2  a13 x3 ...  a1n xn  a11 Из второго уравнения x2: 1 x2  b2  a21 x1  a23 x3 ...  a2 n xn  и так далее. В результате получим систему: a22 © Лекции подготовлены доц. Мусиной М.В. c12 x2  c12 x3  ...  c1n xn  f1 x1   c23 x3  ...  c2 n xn  f2 x2  c21 x1  .  ... xn  c n1 x1  cn 2 x2  cn 3 x3  ...  ...  fn На главной диагонали матрицы С стоят нули, а ненулевые элементы выражают по формулам: a b cij   ij , fi  i ( i, j = 1, 2, … n; i  j) (3) aii aii Для выполнения данного преобразования необходимо чтобы диагональные элементы исходной матрицы не были нулевыми. Если матрица преобразована к виду (2) по данным формулам (3), то метод итераций называют методом Якоби. Для системы (1) метод итерации сходится, если модули диагональных коэффициентов для каждого уравнения системы больше суммы модулей всех остальных коэффициентов (не считая свободных членов). В качестве критерия окончания расчета (достижения заданной точности) можно ис- пользовать простое соотношение X n 1  X n   . Пример. Решим систему уравнений с точностью ε = 0,001 6,25x1  x2  0,5 x3  7,5  - x1  5 x2  2,12 x3  -8,68 0,5x  2,12 x  3,6 x  -0,24 1 2 3  Перепишем систему пересчитав коэффициенты по формулам (3). 0,16 x2  0,08 x3  1,2 x1   0,424 x3  1,786 x2  0,2x1  x  -0,1389x  0,58889 x  0,0667 1 2  3 0.16  0.08    1 .2      Здесь C   0.2  0.424  , F    1.786  .   0.1389  0.5889    0.0667      0    За начальное приближение примем X   0  . 0    Результаты итераций запишем в таблицу. n 0 1 2 3 4 (n) x1 0 1,2 0,9276 0,902 0,8449 (n) x2 0 -1,736 -1,4677 -1,885 -1,8392 x3(n) 0 -0,0667 0,789 0,6688 0,9181 n 1  n  1,736 0,8557 0,4173 0,2493 X X 0  n x1n) x2(n) x3(n) … … … … 12 0,8006 -1,9985 0,9987 13 0,8003 -1,9993 0,9990 14 0,8002 -1,9995 0,9995 15 0,8001 -1,9998 0,9997 © Лекции подготовлены доц. Мусиной М.В. X n 1  X n  … 0,0018 0,0008 0,0005 0,0003 После 15 итераций достигнута искомая точность, и можно записать решение. x1 = 0,8; x2 = -2,0; x3 = 1,0. Методы численной интерполяции и аппроксимации. Процесс научного познания часто приводит к необходимости обработки числовой информации. При обращении с данными научного и инженерного характера особенно важными являются такие вычислительные средства как интерполяция и приближение кривыми (аппроксимация). Некоторые типичные ситуации: 1. Функция f задана таблицей своих значений: yi = f(xi). Требуется найти значение в точке x  xi. 2. Функция f имеет громоздкое выражение и требует для вычисления значительных затрат времени. 3. Значения функции f найдены экспериментально с погрешностью, т.е. yi* вместо yi. Возникающие проблемы решают следующим образом: функцию f(x) заменяют на g(x), значения которой принимают за приближенные значения функции f(x). Линейная интерполяция. Простейшим видом интерполяции является кусочно-линейная интерполяция. Она состоит в том, что заданные точки (xi , yi ), i = 0, ..., n соединяются прямолинейными отрезками, и функция f(x) приближается полученной ломаной. Уравнения каждого отрезка ломаной в общем случае разные. Для каждого из интервалов (xi-1, xi) в качестве уравнения интерполянты используется уравнения прямой, проходящей через две точки. В частности, для i-го интервала можно написать уравнение прямой, проходящей через точки (xi-1 , yi-1 ) и (xi ,yi ), в виде (y - yi-1 )/(yi - yi-1 ) = (x - xi-1 )/(xi - xi-1 ) Отсюда y = aix + bi  g(x), x  (xi-1, xi) ai = (yi -yi-1 )/(xi - xi-1 ), bi = yi-1 - aixi-1 , i = 1, ..., n. Следовательно, при использовании кусочно-линейной интерполяции сначала нужно определить номер i интервала, в который попадает значение аргумента x, затем подставить x и i в формулу и найти приближенное значение функции в точке x. Обычно полагают, что аппроксимируя истинную кривую более сложной линией можно уточнить полученный результат. Например, можно найти многочлен n – ной степени Pn(x), © Лекции подготовлены доц. Мусиной М.В. аппроксимирующий функцию f(x) кривой, проходящей через все n + 1, заданные в таблице точки (xi , yi ), где i = 0, ..., n. В этом случае многочлен должен удовлетворять условиям Pn(xi) = yi, для всех i. Точки (xi , yi ) называют узлами интерполирования, а искомый многочлен – интерполяционным. Интерполяция многочленами. Одна из форм записи интерполяционного многочлена - многочлен Лагранжа. Интерполяционный многочлен для этого метода запишем в виде. Pn(x) = y0b0(x) + y1b1(x) + … + ynbn(x), где все bj(x) – многочлены степени n, коэффициенты которых можно найти с помощью n + 1 уравнений Pn(xi) = yi, i = 0, ..., n. Подставим данные условия в интерполяционный многочлен и получим систему уравнений: y0b0(x0) + y1b1(x0) + … + ynbn(x0) = y0 … y0b0(xn) + y1b1(xn) + … + ynbn(xn) = yn Если значения bj(xi) выбраны так, что 1 при i  j , то уравнения будут удовлетворены. b j x i    0 при i  j Это означает, что любой многочлен bj(x) равен нулю во всех точках xi кроме xj и его можно представить в виде: b j x   C j x  x0 x  x1  x  x j 1 x  x j 1  x  xn  Так как bj(xj) = 1, то коэффициент Cj определяется выражением C j  1 x j  x0 x j  x1  x j  x j 1 x j  x j 1  x j  xn  Т.е многочлен x  x0 x  x1  x  x j 1 x  x j 1  x  xn  bj ( x )  является требуемым многоx j  x0 x j  x1  x j  x j 1 x j  x j 1  x j  xn  членом степени. Подставим все в формулу для Pn(x) и получим. n Pn ( x )   j 0 x  x0 x  x1  x  x j 1 x  x j 1  x  xn  x j y j  Ln x   x0 x j  x1  x j  x j 1 x j  x j 1  x j  x n  Это интерполяционный многочлен Лагранжа. В инженерной практике наиболее часто используется интерполяция многочленами первой, второй и третьей степени (линейная, квадратичная и кубическая интерполяции). При n = 1 – линейная формула. x  x1 x  x0  y1 . (Совпадает с уже полученной формулой линейной инx0  x1 x1  x0 терполяции по двум узлам x0, x1.) L1 x   y 0 При n = 2 L2 x   y 0 x  x1 x  x2   y x  x0 x  x2   y x  x0 x  x1  . x0  x1 x0  x2  1 x1  x0 x1  x2  2 x2  x0 x2  x1  Существует и другой подход для построения интерполяционного многочлена— метод Ньютона (метод разделённых разностей). Но фактически формулы Лагранжа и Ньютона порождают один и тот же полином, разница только в алгоритме его построения. © Лекции подготовлены доц. Мусиной М.В. Аппроксимация функций по методу наименьших квадратов. В инженерной деятельности часто возникает необходимость описать в виде функциональной зависимости связь между величинами, заданными таблично или в виде набора точек (xi , yi ), где i = 0, ..., n. Как правило, эти табличные данные получены экспериментально и имеют погрешности. При аппроксимации желательно получить относительно простую функциональную зависимость (например, многочлен), которая позволила бы «сгладить» экспериментальные погрешности, вычислять значения функции в точках не содержащихся в исходной таблице. Эта функциональная зависимость должна с достаточной точностью соответствовать исходной табличной зависимости. В качестве критерия точности чаще всего используют критерий наименьших квадратов, т.е. определяют такую функциональную зависимость n 2 f(x), при которой R   y i  f ( xi ) обращается в минимум. i 0 Рассмотрим в качестве функциональной зависимости многочлен. n 2 Pm x   a0  a1 x  a2 x 2    am x m , тогда R   y i  Pm ( x i ) . i 0 Условия минимума – нулевые частные производные по всем переменным a0, a1, a2, … n R am. Т.е.   2 y i  a0  a1 xi    am xi xik  0 , или ak i 0 n  y i  a0  a1 x i    am xi xik  0 , k = 0, 1, 2,…, m. i 0 Собираем коэффициенты при неизвестных a0, a1, a2, … am, получаем систему уравнений: n n n n n i 0 i 0 i 0 i 0 i 0 a0  x ik a1  xik 1 a2  x ik 2    am  xik m   xik y i , k = 0, 1, 2,…, m. Можно ввести обозначения: n n i 0 i 0 ck   xik , bk   xik y i и переписать систему в развернутом виде. c0 a0  c1a1  c2 a2   c mam  b0  c1a0  c2 a1  c3 a2   cm 1am  b1  ............................................... cm a0  cm 1a1  cm 2 a2   c2 m am  bm Матрица данной системы называется матрицей Грамма. Решая эту систему линейных уравнений, получаем коэффициенты a0, a1, a2, … am, которые являются искомыми параметрами эмпирической формулы. Рассмотрим два частных случая m = 1 и m = 2. 1. Линейная аппроксимация (m = 1). P1(x) = a0 + a1x © Лекции подготовлены доц. Мусиной М.В. n n n n n i 0 i 0 i 0 i 0 i 0 c0   xi0  n  1 , c1   xi , c2   x i2 , b0   y i , b1   xi y i Таким образом наша система уравнений имеет вид: n n   n  1  a  x a   i 1  yi  i 0 i 0 n n n  x a   x 2 a   x y i 0 i 0 i 0 i 1 i 0 i i Решаем ее методом Крамера. n n n 1  n i 0 i 0  xi  xi2 , a0  n  xi i 0 n i 0 n  xi y i  x  xi i 0 n n  yi i 0  i 0 n n 2 i i 0  yi n 1  xi  xi y i , a1  i 0 i 0  И получаем искомую функцию y = a0 + a1x. 2. Квадратичная аппроксимация (m = 2). P2(x) = a0 + a1x + a2x2 Кроме c0, c1, c2, b0, b1 рассчитываются n n n i 0 i 0 i 0 c3   x i3 , c4   x i4 , b2   y i xi2 .  c0 c1 c2 b0    Расширенная матрица системы уравнений:  c1 c2 c3 b1  , решив которую полуc c c b  3 4 2  2 чим искомые коэффициенты a0, a1, a2. Пример. Получить эмпирическую формулу для функции f(x), заданной таблицей, используя метод наименьших квадратов. x 0,75 1,5 2,25 3 3,75 f(x) 2,3 1,3 1,0 2,2 4,2 Табличные данные изобразим на графике, из которого видно, что в качестве эмпирической функции можно взять параболу. 4,5 4 3,5 3 2,5 2 1,5 1 0,5 0,5 1 1,5 2 2,5 3 3,5 4 © Лекции подготовлены доц. Мусиной М.В. Найдем коэффициенты системы уравнений из таблицы. n x y x2 x3 x4 2,3 0,75 0,5625 0,421875 0,316406 1 1,5 1,3 2,25 3,375 5,0625 2 2,25 1 5,0625 11,39063 25,62891 3 3 2,2 9 27 81 4 3,75 4,2 14,0625 52,73438 197,7539 xy 1,725 1,95 2,25 6,6 15,75 x2y 1,29375 2,925 5,0625 19,8 59,0625 n  i 0 11,25 11 30,9375 94,92188 309,7617 28,275 88,14375 5 a0  11,25 a1  30,94 a2  11  11,25 a0  30,94 a1  94,92a2  28,27 30,94a  94,92 a  309,76 a  88,14 1 2  Решение этой системы a0 = 4,54; a1 = -3,66; a2 = 0,95, и зависимость описывается формулой y = 4,54 – 3,66x + 0,95x2. Погрешность приближения многочленом по методу наименьших квадра1 n y i  Pm xi 2 тов:    n  1 i 0 В случае нашего примера погрешность составит ε = 0,01. Если экспериментальные точки располагаются вдоль некоторой линии, сходной по форме, например, с графиком гиперболической, показательной, логарифмической или других функций с неизвестными параметрами выбирается в качестве аппроксимирующей. Затем проводится линеаризация этой функции с помощью замены переменных и задача сводится к аппроксимации зависимости многочлена первой степени. Например: а). Показательная зависимость: y  ab x , приводится к линейному виду путем логарифмирования ln y  ln a  x ln b . б) Степенная y  ax b , аналогично ln y  ln a  b ln x . x в) Гиперболическая y  , приводится к линейному виду введением новой пеa  bx ременной Y=x/y. Пример. Результаты десяти наблюдений представлены в таблице. Установить вид зависимости между этими величинами и найти параметры эмпирической формулы. x 0 10 11 16 21 27 32 37 43 48 y 8,4 6,2 5,6 5,1 4,2 3,4 3,1 2,5 2,1 1,9 Точечная диаграмма позволяет предположить, что зависимость между x и y экспоненциальная y = aekx. Вводим новые переменные Y =lny, X = x, b = lna, сводим зависимость к линейной Y = kX + b. © Лекции подготовлены доц. Мусиной М.В. 9 8 7 6 5 4 3 2 1 10 20 30 40 50 60 Перестраиваем табличную зависимость. x 10 11 16 21 27 32 37 43 48 lny 2,128 1,825 1,723 1,629 1,435 1,824 1,131 0,916 0,71 0,642 По этим данным методом наименьших квадратов подберем аппроксимирующую линейную функцию. y = a0 + a1x. После решения системы получим a0 = 2,11, a1 = -0,031 = k. Так как a0 =b = lna, то a = ea0 = e2,11 = 8,251. Таким образом, полученная зависимость y = 8,251·e - 0,031·x. Численное интегрирование. К вычислению определенного интеграла сводятся многие практические инженерные и научные задачи. Например, вычисление площадей фигур, объемов тел, работы силы и др. Постановка задачи: Вычислить определенный интеграл b  f ( x )dx , где f(x) – некоторая заданная на отрезке [a, b] функция. a Далеко не всегда можно вычислить интеграл в аналитическом виде по известной из курса математического анализа формуле Ньютона – Лейбница. Иногда даже при известной первообразной ее вид может оказаться очень сложным для вычислений. Также формулу Ньютона – Лейбница нельзя применить, если функция получена экспериментально в виде таблицы. Во всех этих случаях применяют численное интегрирование. На практике наиболее широко используют квадратурные формулы: b n  f ( x )dx   Ai  f xi  i 0 a где xi – точки из отрезка [a, b] (узлы квадратурной формулы), Аi – числовые коэффициенты (веса) квадратурной формулы. n За приближенное значение интеграла принимаем  Ai  f xi  - квадратурную сум- i 0 му. b n R   f ( x )dx   Ai  f xi  - погрешность квадратурной формулы. a i 0 © Лекции подготовлены доц. Мусиной М.В. Для вывода квадратурной формулы отрезок интегрирования [a, b] разбивается на n b равных частичных отрезков [xi - 1, xi], i = 1, …, n длиной h = (b – a) / h, а интеграл  f ( x )dx a заменяется суммой частичных интегралов b n xi  f ( x )dx    f x dx . i  0 xi 1 a Затем подынтегральная функция f(x) на частичном отрезке [xi - 1, xi] заменяется некоторым интерполяционным многочленом невысокой степени m Lm,i(x), и вычисляется интеграл b xi xi xi 1 xi 1  f x dx   Lm,i x dx и приближенное значение интеграла будет получено в виде: n  f ( x )dx   Ai  f xi  . В зависимости от вида интерполяционного многочлена будут поa i 0 лучаться различные квадратурные формулы. Если будет использоваться интерполяционный многочлен нулевой степени, то будет получена формула прямоугольников, если на отрезке [xi - 1, xi] используется линейная интерполяция, то получается формула трапеций, если квадратичная – формула Симпсона. Но можно вывести простейшие квадратурные формулы исходя из геометрического смысла определенного интеграла, а именно находя приближенную площадь криволинейной трапеции. n n i 1 i 1 S   si   f xi xi xi  xi  xi 1 Формула прямоугольников. Можно приближенно заменить площадь криволинейной трапеции площадью ступенчатой фигуры, состоящей из n прямоугольников. Площадь i – того прямоугольника находится как произведение длины отрезка основания [xi - 1, xi] на значение функции в середине отрезка. si  h  f xi 1 2  Суммируя площади всех элементарных трапеций, получим квадратурную формулу средних прямоугольников. © Лекции подготовлены доц. Мусиной М.В. b n     f ( x )dx   f xi 1 2 h h y1 2  y 3 2  ...  y n 3 2  y n 1 2 a  i 1  si  h  f xi   h  f xi 1   h  f xi 1 2  Иногда используют формулы правых и левых прямоугольников: b n  f ( x )dx   f xi h hy1  y 2  ...  y n 1  y n  a b i 1 n  f ( x )dx   f xi 1 h hy 0  y1  ...  y n 2  y n 1  a i 1 Геометрическая иллюстрация этих формул. Формула трапеций. Формулу трапеций можно также получить из геометрических соображений. Отрезок ba интегрирования [a, b] разбивается на n частей длины h  . ( h - шаг разбиения.) n x1 , x2 ,  x n - абсциссы точек деления и соответствующие ординаты кривой y1, y2, …, yn соединяются ломаной линией. В результате построения криволинейная трапеция разбивается на ряд вертикальных полосок одной и той же ширины h, каждую из которых приближенно h можно принять за трапецию. Площадь элементарной трапеции: si  f xi 1   f xi  2 Суммируя площади этих трапеций, получаем формулу: b y  y  f ( x )dx  h 20  y1  y 2    y n1  2n  (Формула трапеций.) a Формула Симпсона. Более точную формулу можно получить, если профиль криволинейной полоски считать параболой (используя многочлен Лагранжа второй степени), а не прямой линией как в формуле трапеций. В этом случае можно получить формулу Симпсона: © Лекции подготовлены доц. Мусиной М.В. b  f ( x )dx  a ba  f ( a )  2 f (a  2h )    f (a  (n  2)h )  3n  4 f ( a  h )    f (a  ( n  1)h   f (b), h  ba . n Для практической оценки погрешности интеграла используется правило Рунге. Для этого проводят вычисления с шагом h и h/2 и получают приближенные значения интегралов Ih и Ih/2. За погрешность приближенного значения интеграла принимают величину: Ih 2  Ih / 3 - для формулы прямоугольников и трапеций, Ih 2  Ih / 15 - по формуле     Симпсона. В адаптивных алгоритмах автоматически определяют величину шага h, таким образом, чтобы результат удовлетворял заданной точности. За окончательные значения интеграла принимают: Ih 2  Ih 2  Ih / 3 - по формуле прямоугольников, Ih 2  Ih 2  Ih / 3 -       формуле трапеций, Ih 2  Ih 2  Ih / 15 - формуле Симпсона. Решение обыкновенных дифференциальных уравнений. Инженеру часто приходится иметь дело с техническими системами и технологическими процессами, характеристики которых непрерывно меняются со временем t. Эти явления, как правило, подчиняются физическим законам, которые формулируются в виде дифференциальных уравнений. Например, кинетика химических реакций, динамика биологических популяций, движение космических объектов, модели экономического развития, все эти процессы исследуются с помощью обыкновенных дифференциальных уравнений (ОДУ). Обыкновенным дифференциальным уравнением первого порядка называется уравнение вида F(x,y,у')=0 или у'=f(x,y). Функция y(x), при подстановке которой уравнение обращается в тождество, называется решением дифференциального уравнения. Для решения обыкновенных дифференциальных уравнений могут применяться точные методы решения, а также приближенные (аналитические) методы, такие, например, как разложение в числовой ряд, рассмотренные в курсе высшей математики. К приближенным методам решения относятся и численные методы. При решении ДУ численным методом приближенные значения искомой функции y(x) находят для ряда значений x из отрезка, на котором нужно найти решение. С помощью численных методов невозможно найти общее решение, а только частное решение, то есть решение задачи Коши. Подавляющее большинство возникающих на практике задач Коши невозможно решить без использования вычислительной техники. Поэтому в инженерных и научно-технических расчетах численные методы решения задачи Коши играют важную роль. Рассмотрим несколько численных методов решения дифференциальных уравнений первого порядка. Постановка задачи: Решить уравнение y' = f(x, y), при начальном условии y(x0) = y0, где x, x0[a, b]. Первый этап на пути построения численного методов решение на отрезке [a, b] задачи Коши состоит в замене отрезка [a, b] – области непрерывного изменения аргумента x – множеством конечного числа точек a=x0 < x1 < ... < xi < ...< xn=b, называемого сеткой. Сами точки xi - узлы сетки. Величина h =xi - xi-1 - шаг сетки. Будем рассматривать сетки, в © Лекции подготовлены доц. Мусиной М.В. которых шаг h - постоянный. В таком случае xi = a+ i h, h=(b-a)/ n, и сетка называется равномерной. Следующий этап в построении численного метода состоит в замене задачи Коши y' = f(x, y), y(x0) = y0 ее дискретным аналогом – решая которую можно последовательно найти таблицу приближенных значений y0, y1, ..., yi, ... yn, решения y(x) в узлах сетки, то есть y(xi)= yi. Численный метод решения задачи Коши называется одношаговым, если для вычисления решения в точке xi + h используется информация о решении только в точке xi. Простейший одношаговый метод численного решения задачи Коши – метод Эйлера. Метод Эйлера. Метод Эйлера является сравнительно грубым и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других методов. Идея метода заключается в том, что на малом промежутке изменения независимой переменной интегральная кривая дифференциального уравнения заменяется отрезком прямой (касательной.) Пусть необходимо найти решение уравнения y' = f(x, y) (1) с начальным условием y(x0) = y0 (2). Разложим искомую функцию y(x) в ряд вблизи точки x0 и ограничимся первыми двумя членами разложения y(x) = y(x0) + y´(x0)(x –x0) + …. Учитываем уравнение (1) и обозначив h = x –x0, получаем y(x) ≈ y(x0) + f(x0, y0)·h. Эту формулу можно применять многократно, находя значения функции во все новых и новых точках. уi+1 = уi+h·f(xi, yi) (i=0,1,2…). (3) Такой метод решения обыкновенных дифференциальных уравнений называется методом Эйлера. Запускается метод из начальных условий y(x0) = y0. Геометрически метод Эйлера означает, что на каждом шаге мы аппроксимируем решение (интегральную кривую) отрезком касательной, проведенной к графику решения в начале интервала. При этом искомая интегральная кривая у=у(х), проходящая через точку М0(х0, у0), заменяется ломаной М0М1М2… с вершинами Мi(xi, yi) (i=0,1,2,…); каждое звено МiMi+1 этой ломаной, называемой ломаной Эйлера, имеет направление, совпадающее с направлением той интегральной кривой уравнения, которая проходит через точку Мi. Метод Эйлера имеет низкую точность и для того чтобы получить приемлемую погрешность необходимо делать шаг расчета достаточно маленьким. Можно улучшить точность метода Эйлера, улучшив аппроксимацию производной. Это можно сделать, например, используя среднее значение производной в начале и в конце интервала. Модифицированный метод Эйлера. Рассматриваем дифференциальное уравнение y' = f(x, y) с начальным условием y(x0) = y0. Вычислим значение функции в следующей точке по методу Эйлера: y*i+1 = уi+h·f(xi, yi) и используем это значение для нахождения производной в конце интервала f(xi+1, y*i+1). Вместо значения f(xi, yi) в формуле Эйлера используем полусумму f(xi, yi) и f(xi+1, y*i+1). уi+1 = уi+h·½·(f(xi, yi) + f(xi+1, y*i+1)). (4) © Лекции подготовлены доц. Мусиной М.В. Разложим функцию в ряд Тейлора вблизи точки x0 и добавим еще один член разложения: y(x) = y(x0) + y´(x0)(x –x0) +½ y´´(x0)(x –x0)2 + …. Заменим вторую производную разy x0  h   y x0  ностным уравнением y ( x0 )  и подставим в разложение. h y(x + h) ≈ y(x0) + y´(x0) ·h +½ ( y´(x0 + h) - y´(x0) )/ h ·h 2 Преобразуем y(x + h) ≈ y(x0) + y´(x0) ·h +½ ( y´(x0 + h) - y´(x0) ) ·h y(x + h) ≈ y(x0) + ½ ( y´(x0 + h) + y´(x0) ) ·h Эта формула совпадает с уже полученным выражением (4). При этом точное значение y´(x0 + h) заменено приближенным f(xi+1, y*i+1). При вычислении модифицированным методом Эйлера, сначала вычисляют значение y*i+1, а затем уже подставляют его в формулу (4). Метод Рунге – Кутта. Методом Рунге – Кутта четвертого порядка точности называют одношаговый метод, относящийся к широкому классу методов Рунге-Кутты. В этом методе величины yi+1 вычисляются по следующим формулам: h y i 1  y k  m1  2 m2  2 m3  m4   , 6 где m1  f xi , y i  h h  m2  f  xi  , y i  m1   2 2  h h  m3  f  xi  , y i  m2   2 2  h   m4  f  xi  , y i  m3  h  . 2   Пример. Решить дифференциальное уравнение y' = y + 3x, при начальном условии y(0) = -1, где x[0; 0,5]. Точное аналитическое решение задачи Коши: y = 2ex – 3x – 3. Решим методом Эйлера. Выберем для решения шаг h = 0,1. Формула для расчета уi+1 = уi+h·f(xi, yi), т.е. уi+1 = уi+h·(3xi + yi), xi = 0,1· i Чтобы получить решение в точке x1, используем значения в начальной точке. у1 = у0+h·(3x0 + y0), у1 = - 1+0,1·(3·0 - 1) = - 1,1, x1 = 0,1· 1 = 0,1. Следующая точка: x2 = 0,1· 2 = 0,2 y2 = у1+h·(3x1 + y1), у1 = - 1,1+0,1·(3·0,1 - 1,1) = - 1,18 Запишем результаты расчетов в таблицу. i x 1 2 3 4 5 0,1 0,2 0,3 0,4 0,5 Точное решение -1 -1,08966 -1,15719 -1,20028 -1,21635 -1,20256 h =0.1 -1 -1,1 -1,18 -1,238 -1,2718 -1,27898 Относительная погрешность 0,940167 1,932671 3,046657 4,359915 5,975273 © Лекции подготовлены доц. Мусиной М.В. Проведем компьютерный расчет с шагами h = 0,01 и h = 0,001. Результаты показаны в таблице. x 0,1 0,2 0,3 0,4 0,5 Точное решение -1 -1,08966 -1,15719 -1,20028 -1,21635 -1,20256 h=0,01 -1 -1,09076 -1,15962 -1,2043 -1,22227 -1,21074 Относительная погрешность 0,100728 0,209596 0,334903 0,48686 0,680125 h=0,001 -1 -1,08977 -1,15744 -1,20069 -1,21695 -1,20338 Относительная погрешность 0,010135 0,021095 0,033714 0,049022 0,068496 Оценка погрешности по правилу Рунге. Численное решение дифференциального уравнения всегда получается с погрешностью. Практически оценить погрешность метода позволяет правило Рунге. Сначала вычисляют приближенное решение с шагом h, затем - с шагом h/2. и сравнивают величины y ih  max y ih  y ih 2  h 2 и yi .Оценкой погрешности метода будет величина: , где p – порядок 2p 1 точности метода. Для метода Эйлера (p = 1): max y ih  y ih 2 Для модифицированного метода Эйлера(p = 2): Для метода Рунге-Кутта (p = 4): h i max y  y 15 h2 i . max y ih  y ih 2 3