Математика. Численные методы

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ «ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ – УЧЕБНО-НАУЧНОПРОИЗВОДСТВЕННЫЙ КОМПЛЕКС» Г.А. Семенова, Т.В. Савостикова МАТЕМАТИКА ЧИСЛЕННЫЕ МЕТОДЫ Рекомендовано ФГБОУ ВПО «Госуниверситет-УНПК» для использования в учебном процессе в качестве учебного пособия для высшего профессионального образования Орёл 2015 1 УДК 519.2 (075) ББК 22.17я7 Рецензенты: доктор технических наук, профессор, заведующий кафедрой «Высшая математика» Федерального государственного бюджетного образовательного учреждения высшего профессионального образования «Государственный университет – учебно-научно-производственный комплекс» Владимир Александрович Гордон, кандидат педагогических наук, доцент кафедры «Алгебра и математические методы в экономике» Федерального государственного бюджетного образовательного учреждения высшего профессионального образования «Орловский государственный университет» Ирина Ивановна Чернобровкина Семенова Г.А. С-17 Математика. Численные методы: учебное пособие для высшего профессионального образования / Г.А. Семенова, Т.В. Савостикова. – Орёл: ФГБОУ ВПО «Госуниверситет-УНПК», 2015. – 56 с. Учебное пособие содержит теоретический материал и примеры решения задач по основам численных методов. Предназначено студентам укрупнённой группировки направлений «Информатика и вычислительная техника» подготовки бакалавров очной формы обучения, изучающим дисциплину «Математика». УДК 519.2 (075) ББК 22.17я7 © ФГБОУ ВПО «Госуниверситет-УНПК», 2015 2 СОДЕРЖАНИЕ ВВЕДЕНИЕ ................................................................................................. 4 1. Приближённые числа и действия над ними ........................................ 5 1.1 Погрешность ...................................................................................... 5 1.2 Абсолютная и относительная погрешность ................................... 6 1.3 Десятичная запись приближенных чисел ....................................... 7 1.4. Действия над приближенными числами ....................................... 8 1.5. Строгий учёт погрешности по методу (способу) границ .......... 11 2. Классификация численных методов .................................................. 12 3. Аппроксимация функции .................................................................... 14 3.1. Постановка задачи.......................................................................... 14 3.2. Интерполяционный многочлен Лагранжа. .................................. 15 3.3. Интерполяционный многочлен Ньютона для равноотстоящих узлов .................................................................. 20 3.4. Аппроксимация функции. Метод наименьших квадратов ........ 26 3.5. Нелинейная аппроксимация .......................................................... 31 4.Решение нелинейных уравнений ......................................................... 33 4.1 Постановка задачи........................................................................... 33 1.2. Основные этапы отыскания решения .......................................... 34 1.3. Метод половинного деления (метод бисекции).......................... 34 4.4 Метод касательных Ньютона ......................................................... 35 4.5. Метод секущих (метод хорд) ........................................................ 37 5. Численное интегрирование функции одной переменной ............... 41 5.1 Постановка задачи........................................................................... 41 5.2 Метод прямоугольников ................................................................ 41 5.3. Метод трапеций .............................................................................. 43 5.4. Метод Симпсона (метод парабол) ................................................ 44 5.5. Правило Рунге практической оценки погрешности ................... 46 6. Численное решение обыкновенных дифференциальных уравнений .................................... 48 6.1. Постановка задачи.......................................................................... 48 6.2. Метод Эйлера ................................................................................. 50 6.3. Метод Рунге-Кутта ......................................................................... 51 ЛИТЕРАТУРА .......................................................................................... 56 3 ВВЕДЕНИЕ Исследования различных явлений или процессов математическими методами осуществляются с помощью математических моделей. Математическая модель представляет собой формализованное описание на языке математики исследуемого объекта. Это может быть система линейных, нелинейных или дифференциальных уравнений, система неравенств, определённый интеграл, многочлен с неизвестными коэффициентами и т.д. Построенная модель должна охватывать важнейшие характеристики исследуемого объекта и отражать связи между ними. После того, как математическая модель составлена, переходят к постановке вычислительной задачи. При этом устанавливают, какие характеристики математической модели являются исходными (входящими) данными, какие – параметрами модели, а какие – выходными данными. Проводится анализ полученной задачи с точки зрения существования и единственности решения. На следующем этапе выбирается метод решения задачи. Во многих конкретных случаях найти решения задачи в явном виде не представляется возможным, так как оно не выражается через элементарные функции. Такие задачи можно решить лишь приближённо. Под вычислительными (численными) методами подразумеваются приближенные процедуры, позволяющие получать решения в виде конкретных числовых значений. Вычислительные методы, как правило, реализуется на ЭВМ. Для решения одной и той же задачи могут быть использованы различные вычислительные методы, поэтому нужно уметь оценивать качество различных методов и эффективность их применения для данной задачи. Затем для реализации выбранного вычислительного метода составляется алгоритм и программа на ЭВМ. Результаты расчёта анализируются и интерпретируются. При необходимости корректируются параметры метода, а иногда математическая модель, и начинается новый цикл решения задачи. 4 1. Приближённые числа и действия над ними 1.1 Погрешность Численное решение любой задачи осуществляется приближенно, с различной точностью. Главная задача численных методов – фактическое нахождение решения с требуемой или, по крайней мере, оцениваемой точностью. Отклонение истинного решения от приближённого называется погрешностью. Полная погрешность вычисления состоит из двух составляющих: неустранимая и устранимая погрешности. Неустранимая погрешность обусловлена неточностью исходных данных и никоим образом не может быть уменьшена в процессе вычислений. Устранимая погрешность в свою очередь состоит из погрешности аппроксимации (метода) и погрешности вычислений. Существуют четыре источника погрешностей, возникающих в результате численного решения задачи. 1. Математическая модель. Погрешность математической модели связана с её приближенным описанием реального объекта. Погрешность математической модели неустранима. 2. Исходные данные. Исходные данные обычно содержат погрешности, так как они либо неточно измерены, либо являются решением некоторых вспомогательных задач. В среднем погрешность измерений составляет 1-10%. Погрешность исходных данных считается неустранимой. 3. Метод вычислений. Применяемые для решения задачи методы, как правило, являются приближёнными. Погрешность метода можно оценить и проконтролировать. Следует выбирать погрешность метода так, чтобы она была не более чем на порядок меньше неустранимой погрешности. 4. Округление в вычислениях. Погрешность округления возникает из-за того, что вычисления проводятся с конечным числом значащих цифр. При решении большѝх задач производятся миллиарды вычислений, но так как погрешности имеют разные знаки, то они частично взаимно компенсируются. 5 1.2 Абсолютная и относительная погрешность Пусть А – точное значение некоторой величины, а – её приближенное значение. Пишут А ≈ а. Например, A  2  1,41  a . Если а > A, то а называется приближением по избытку, если же а < A – приближением по недостатку. Абсолютной погрешностью приближения а называется число ∆ = |A – a|. Его достаточно знать для характеристики отличия приближённого значения от точного. Однако, часто такое точное значение исследуемой величины неизвестно, а значит, неизвестно и ∆. Поэтому целесообразно ввести оценку абсолютной погрешности. Предельной абсолютной погрешностью  a называется точная верхняя грань абсолютной погрешности ∆, то есть  a  sup{} . Отсюда,   A  a  a  a  a  A  a  a  A  a  a . При этом a   a оказывается приближением по недостатку, a   a – приближением по избытку. Пример. Рассмотрим число A    3,14  a . Известно, что 3,14    3,15 . Следовательно, абсолютная погрешность нашего приближения     a  0,01, в качестве предельной абсолютной погрешности можно взять число  a  0,01. Если же учесть, что число   (3,140; 3,142) , то получим лучшую оценку  a  0,002 , а значит, можно записать   3,140  0,002 . Для характеристики качества измерения величины А одной абсолютной погрешности мало. Например, если измерить длину дачного участка и расстояние от Орла до Москвы с погрешностью 0,1 м, то ясно, что точность второго измерения будет выше.  Относительная погрешность приближения а – это число   . a Предельной относительной погрешностью приближения а называется положительное число  a такое, что   a . Очевидно, что за предельную относительную погрешность можно принимать    а  a  a 100% . a a 6 0,002 100%  3,140  0,00064 100%  0,064%  0,01% – Приемлемая для инженеров точность. Упражнение. Найдите предельную относительную погрешность следующих приближений A1  34,8  0,2 ; A2  0,00464  0,00004 ; A3  4327  30 . Пример. A    3,140  a ,  a  0,002 =>  а  1.3 Десятичная запись приближенных чисел Форма записи A  a  a приближенных чисел не всегда удобна. Целесообразно иметь такую форму, способ записи которой бы позволял по десятичной иметь такую форму, способ записи которой бы позволял по десятичной записи числа установить абсолютную погрешность приближения. Цифра α в десятичной записи приближённого числа называется верной в широком смысле, если абсолютная погрешность приближения не превосходит единицы того разряда, которому принадлежит цифра α. Если же абсолютная погрешность приближения a не превосходит 1 2 единицы указанного разряда, то её (α) называют верной в узком смысле. Пример. Рассмотрим число A  7,158  0,0009 . Абсолютная погрешность  a  0,0009  0,001 (три цифры после запятой), следовательно, все цифры приближения a  7,158 верны в широком смысле. 1 С другой стороны,  a  0,0009   0,001  0,0005 , поэтому α = 8 (тре2 тья цифра после занятой) не является верной в узком смысле. Значащими цифрами называются все цифры числа, начиная от первой слева ненулевой цифры. Например, у числа 0,037 две значащие цифры, у числа 14,80 – четыре. Правила округления, обеспечивающие все верные знаки в узком смысле. 1. Если первая из отбрасываемых цифр меньше пяти, то сохраняемые цифры остаются без изменения. 2. Если первая из отбрасываемых цифр равна пяти, а среди следующих за ней цифр есть отличные от нуля, то последнюю из сохраняемых цифр увеличивают на единицу. 7 3. Если первая из отбрасываемых цифр равна пяти, а все следующие за ней являются нулями, то последнюю из сохраняемых цифр увеличивают на единицу, когда она нечетная, и сохраняют неизменной, когда она чётная. 4. Если первая из отбрасываемых цифр больше пяти, то последнюю из сохраняемых цифр увеличивают на единицу. Пример.   3,14159265...  3,1416 . Все пять цифр этого приближения верны и в узком, и в широком смысле. При записи приближённых чисел следует писать только верные значащие цифры, если погрешность не задана другим способом. Изменения формы записи числа, например, на число с плавающей запятой, не должно менять количество верных цифр. Любое число d в десятичной системе можно записать в форме с плавающей запятой, то есть в виде d   m 10 p , где m – мантисса, p – порядок. Способов такой записи много. Среди них выделяют нормальную и стандартную формы записи. Если m < 1 и первая цифра в m после запятой отлична от нуля, то это нормальная форма записи числа. Если же m  1; 10 , то форма записи называется стандартной. Например, 123,4  0,1234 103 – нормальная 123,4  1,234 102 – стандартная форма. форма записи, 1.4. Действия над приближенными числами Теорема. Пусть функция f ( x1 ,.., xn ) дифференцируема в области G  n , a её переменные вычислены с предельными абсолютными погрешностями  x k , k = 1,…,n. Тогда для предельной абсолютной погрешности  f значения функции f справедливо равенство f   xk ,  x k k 1 а для предельной относительной погрешности  f значения функции f справедливо равенство n f  n f   k 1 1 f    xk . f xk 8 Здесь функция f и её частные производные f вычислены в точке xk M ( x1* ,..., xn * ) . Пример. Период колебаний математического маятника вычисляется по формуле T  2  l g . Найдём погрешность вычисления периода, если   3,142 , l  1,2000 м, g  9,8132 м/с2. Решение. В соответствии с договорённостями все знаки считаются верными в узком смысле, то есть    0,00041, l  0,0005 ,  g  0,00005 . По предыдущей теореме l  l T  2     l  2 3  (3 / 2)   g  g g l g 1,2000 3,142  0,0005 1,2000  0,00041   2  3,142   0,00005 ; 3 9,8132 1,2000  9,8132 9,8132 T  0,00028675  0,0004578  0,00000559  0,0003  0,0005 ; 0,0001  T  0,0005 . Следовательно, вычислить значение периода Т колебаний в данных условиях можно только с тремя верными в узком смысле знака1,2000 ми, то есть T  2  3,142   2,197 (1/с). 9,8132 Таким образом, T  2,197  0,001. T  2 Следствие. Предельная абсолютная погрешность суммы равна сумме предельных абсолютных погрешностей слагаемых, то есть n  a1 ... an    ak . k 1 Следствие. Предельная абсолютная погрешность произведения двух функций удовлетворяет соотношению  u v   u  v   v  u . Следствие. Предельная относительная погрешность произведения функций равна сумме предельных относительных погрешностей множителей, то есть 9 n  a1 a2 ...an    ak . n 1 Следствие. Предельная относительная погрешность частного двух функций равна сумма предельных относительных погрешностей, то есть  u v  u  v . Следствие. Предельная относительная погрешность степени  a  в  раз больше предельной относительной погрешности основания, то есть  a      В частности, извлечение корня n-ой степени уменьшает предельную относительную погрешность в n раз. Пример. Даны три числа a = 1,23; b = 3,45; c = 0,15, все знаки которых верны в широком смысле, то есть  a  0,01; b  0,01; c  0,01. а) Найдём сумму чисел S  a  b  c . Решение. Приближённое значение s  1,23  3,45  0,15  4,83 , его абсолютная погрешность  s   a  b  c   s  0,01  0,01  0,01  0,03  S  s   s  4,83  0,03  S [4,80; 4,86]  S  4,8  0,06 . б) Найдём произведение чисел P  a  b . Решение. Приближённое значение p  1,23  3,45  4,2435 . Относительная погрешность множителей 0,01 0,01 a   0,00813; b   0,00290 1,23 3,45 следовательно, относительная погрешность произведения ab  a  b  0,00813  0,01103   ab  p  ab  0,0468  0,05 , то есть у произведения всего один верный в узком смысле знак после запятой. Точное значение P  a  b   ab ; P  4,2435  0,04682  P  4,1967; 4,2903. С учётом вывода о количестве верных знаков P  4,2  0,1. в) Найдём частное чисел D  b c . Решение. Приближённое значение d  3,45 0,15  23 , его относительная погрешность 10  bc  следовательно,  b c  d   b c 0,01 0,01   0,0696 , 3,45 0,15  23  0,0696  1,61  2 . Отсюда, D  23  1,61  D  21,39; 24,61. С учётом вывода о количестве верных знаков D  23  2 . Упражнение. Оценить сверху относительную погрешность при 1 вычислении площади треугольника по формуле S  ab  sin  , в ко2 тором a = 27,3; b = 44,6;   6426  2 . 1.5. Строгий учёт погрешности по методу (способу) границ Если про некоторое число Х известно, что X  a; b, то а – нижняя, b – верхняя граница числа Х. Обозначим нижнюю границу числа Х – НГХ, а верхнюю границу – ВГХ. Вычисление со строгим учетом погрешностей по методу границ позволяют дать точную оценку погрешности полученного результата. В основе метода лежат ранее сформулированные следствия и теоремы, описанные ниже. Теорема 1 НГ X Y  НГ X  НГ Y , ВГ X Y  ВГ x  ВГ y . Теорема 2 НГ X Y  НГ X  НГ Y , ВГ X  y  ВГ x  ВГ y . Теорема 3 Если x, y>0, то НГ XY  НГ X  НГ Y , ВГ XY  ВГ x  ВГ y . Теорема 4 НГ X ВГ X Если x, y>0, то НГ X Y  , ВГ X Y  . НГ Y ВГ Y Замечание. На основании свойств неравенств очевидно, что нижнюю границу можно и нужно округлять по недостатку, а верхнюю – по избытку. Пример. Найдём верхнюю и нижнюю границы для числа B D  C , где B  b  3,45; C  c  0,15 . Решение. Абсолютные погрешности чисел В и С  B  0,01; C  0,01 , следовательно, B [3,44; 3,46], C [0,14; 0,16] . Отсюда 11 НГ В 3,44 ВГ B 3,46   21,5; ВГ B C    24,72 ВГ С 0,16 НГ C 0,14  D [21,5; 24,72]. НГ B C  Замечание. Если для числа Х известны верхняя и нижняя границы, то за приближённое значение величины Х удобно брать их среднее арифметическое ВГ X  НГ X , 2 a в качестве предельной абсолютной погрешности – x ВГ X  НГ X . 2 В рассмотренном выше примере b 21,5  24,72 24,72  21,5 d   23,11,  B C   1,61  D  23,11  1,61 . c 2 2 x  2. Классификация численных методов Под численными методами понимаются методы, которые используются в вычислительной математике для преобразования задач к виду, удобному для реализации на ЭВМ. Различают два класса методов. Прямые методы. Это методы, которые позволяют получить решение после выполнения конечного числа элементарных операций. Иногда прямые методы называют точными, так как при отсутствии ошибок в исходных данных и при выполнении элементарных операции, результат будет точным. Однако, при реализации метода на ЭВМ неизбежны ошибки округления и, как следствие, наличие вычислительной погрешности. Итерационные методы. Суть методов состоит в построении последовательных приближений к решению задачи. Суть методов состоит в построении последовательных приближений к решению задачи. Вначале выбирают одно или несколько начальных приближений, а затем последовательно, используя найденные ранее приближения и однотипную процедуру расчёта, строят новые приближения. В результате такого итерационного процесса можно теоретически построить бесконечную последовательность приближений к решению. Если эта последовательность сходится (что бывает не всегда), то го12 ворят, что данный метод сходится. Отдельный шаг итерационного процесса называется итерацией. На практике вычисления не могут продолжаться бесконечно долго. Поэтому необходимо выбирать критерий окончания итерационного процесса. Критерий окончания связан с требуемой точностью вычислений, то есть вычисления заканчиваются тогда, когда погрешность приближения не превышает заданной величины. Оценки погрешности приближения, полученные до вычислений, называются априорными оценками, а соответствующие оценки, полученные в ходе вычислений, называются апостериорными оценками. Важной характеристикой итерационных методов является скорость сходимости метода. Говорят, что метод имеет р-ый порядок сходимости, если | xn 1  x* | C | xn  x* | p , где xn , xn1 – последовательные приближения, полученные в ходе итерационного процесса вычислений, x* – точное решение, C – константа, не зависящая от n. Говорят, что метод сходится со скоростью геометрической прогрессии со знаменателем q<1, если для всех n справедливо неравенство | xn  x* | Cq n . Итерационный процесс называется одношаговым, если для вычисления очередного приближения xn используется только одно предыдущее приближение xn1 . Итерационный процесс называется kшаговым, если для вычисления xn используется k предыдущих приближений xn 1,..., xn  k . . К решению задач приближенными численными методами предъявляются требования корректности метода. Это достаточно естественные требования существования, единственности и устойчивости решения. В частности, решение задачи называется устойчивым, если оно зависит от исходных данных непрерывным образом, то есть малому изменению условий задачи соответствует малое изменение решения задачи. Если хотя бы одно из условий не выполнено, задача называется некорректной, а значит, необходимо пересмотреть построенную математическую модель. 13 3. Аппроксимация функции 3.1. Постановка задачи Задача приближения (аппроксимации) функции f ( x)  A заключается в том, чтобы для данной функции построить другую, отличную от неё функцию ( x)  B , значение которой достаточно близки к значениям данной функции. Наиболее типичные случаи, приводящие к необходимости решения такой задачи: 1) функция задана таблично в конечном множестве точек, а вычисления нужно произвести в других точках; 2) функция задана аналитически, но её вычисление по формуле затруднительно. Функция f (x) из класса A называется аппроксимируемой, а функция (x) из класса B называется аппроксимирующей (это чаще всего многочлен или рациональная функция). Теорема Вейерштрасса (об аппроксимации). Если функция f (x) непрерывна на отрезке [a,b], то (  0) (( x)  Pn ( x)) ( x [a, b]) [ | f ( x)  Pn ( x) | ], где Pn (x) – многочлен n-ой степени, то есть у любой непрерывной функции есть “близкий” ей многочлен. Замечание. Теорема Вейерштрасса говорит о существовании многочлена, но ничего не говорит о его построении и оценке приближения. Если в качестве аппроксимирующей функции (x) выбрать функцию, значения которой (n  1) -ой точке a  x0  x1  ...  xn 1  xn  b совпадает со значениями функции f (x) , то есть для  i  0,..., n ( xi )  f ( xi ) , то такая аппроксимация называется интерполяцией. Точки x0 ,..., xn называются узлами интерполяции. Разумно ожидать выполнения приближенного равенства ( x)  f ( x)  x [a, b]. Если приближение выходит за интервал [a,b], то речь идёт об экстраполяции функции. Если в качестве критерия близости исходной f (x) и приближенной (x) функций выбрать минимизацию среднеквадратичного отклонения, то получим метод наименьших квадратов. 14 3.2. Интерполяционный многочлен Лагранжа Пусть функция y  f (x) определена на отрезке [a,b] (рис. 1), и известны значения этой функции в узлах a  x0  x1  ...  xn  b . Обозначим эти значения следующим образом: yi  f ( xi ), i  0,..., n . Требуется найти такой многочлен Lm ( x)  a0  a1x  a2 x 2  ...  am x m , который в узлах интерполяции принимает те же значения, что и исходная функция y  f (x) , то есть (1) Lm ( xi )  yi , i = 0,…,n. Многочлен Lm (x), удовлетворяющий описанному условию, называется интерполяционным многочленом. Рисунок 1. Задача интерполяции Другими словами, ставится задача построения функции y  Lm (x), график которой проходит через заданные точки ( xi , yi ), i  0,..., n . Запишем равенства (1) в виде системы (п+1)-го уравнения с (m+1)-ой неизвестной с коэффициентами многочлена a0 , a1,..., am : a0  a1 x0  a2 x02  ...  am x0m  y0 ;  (1’) ...  2 m a0  a1 xn  a2 xn  ...  am xn  yn . Для существования единственного решения необходимо, чтобы число уравнений совпадало с числом неизвестных, то есть m = n. 15 Кроме того, в силу теоремы Крамера, определитель системы (1’) должен быть отличен от нуля. Действительно, определитель системы 1 x0 x02 ... x0n W  ... ... ... ... ... – 1 xn xn2 ... xnn это определитель Вандермонда. Так как по построению xi  x j ;  i, j , то W  0. Итак, система (1’) при n = m имеет единственное решение, что в свою очередь говорит о корректности поставленной задачи. Имеются различные формы записи интерполяционного многочлена. Лагранж предложил записать многочлен Ln (x) в виде линейной комбинации многочленов той же степени n, то есть 1, i  k ; (2) Ln ( x)   yi  qi ( x), где qi ( xk )   ik   0, i  k . i 0 Такое дополнительное условие обеспечивает выполнение равенств Ln ( xi )  yi . В качестве функций qi (x) можно взять функции вида ( x  x0 )...( x  xi 1 )( x  xi 1 )...( x  xn ) , qi ( x)  ( x  x0 )...( x  xi 1 )( x  xi 1 )...( x  xn ) тогда интерполяционный многочлен Лагранжа имеет вид n ( x  x0 )...( x  xi 1 )( x  xi 1 )...( x  xn ) . (3) Ln ( x)   yi  ( x  x )...( x  x )( x  x )...( x  x ) i 0 i 1 i 1 n n Частные случаи интерполяционного многочлена Лагранжа. 1. При n=1 получаем случай линейной интерполяции. x  x0 x  x1 L1 ( x)  y0  y1 (4) x0  x1 x1  x0 График аппроксимирующей функции y  L1 ( x) (рис. 2) представляет собой прямую, проходящую через точки ( x0 , y0 ) и ( x1 , y1 ) . 2. При n=2 получаем случай квадратичной интерполяции. L2 ( x)  ( x  x0 )( x  x2 ) ( x  x0 )( x  x1 ) ( x  x1 )( x  x2 ) y0  y1  y 2 (5) ( x0  x1 )( x0  x2 ) ( x1  x0 )( x1  x2 ) ( x2  x0 )( x2  x1 ) 16 Рисунок 2. График линейной интерполяции Рисунок 3. График квадратичной интерполяции График аппроксимирующей функции y  L2 ( x) (рис. 3) представляет собой параболу, ось симметрии которой перпендикулярна оси Оx, и которая проходит через точки ( x0 , y0 ) , ( x1 , y1 ) , ( x2 , y 2 ) . Пример. Экспериментально получены четыре точки xk yk 23,30 299 24,25 25,25 26,10 328 373 415 Найдём значение функции y  f (25) , используя линейную и квадратичную интерполяцию. Решение. 25  26,10 25  25,25 y (25)  L1 (25)   3,73   415  360,6470588  25,25  26,10 26,10  25,25  360,6471. 17 (25  25,25)( 25  26,10)  328  (24,25  25,25)( 24,25  26,10) (25  24,25)( 25  26,10)   373  (25,25  24,25)( 25,25  26,10) (25  24,25)( 25  25,25)   415  361,3028167  361,3029 . (26,10  24,25)( 26,10  25,25) y (25)  L2 (25)  Следует отметить, что интерполяционный многочлен (3) малопригоден для быстрых вычислений при n>3. На практике используется расчётная схема. Расчётная (практическая) схема вычисления полинома Лагранжа. Запишем интерполяционный многочлен Лагранжа в более удобном виде Ln ( x)  A0  A1 ( x  x0 )  A2 ( x  x0 )( x  x1 )  ... (6)  An ( x  x0 )( x  x1 )...( x  xn 1 ). Тогда Ln ( x0 )  A0  y0 . Рассмотрим теперь другой многочлен L ( x)  y 0 B( x)  n  A1  A2 ( x  x1 )  ...  An ( x  x1 )...( x  xn1 ) . x  x0 y  y0 Тогда B( x1 )  A1  1 . x1  x0 Таблица 1. Расчётная схема вычисления многочлена Лагранжа xk x0 yk y0 x1 y1 x2 y2 x3 y3 … … xn yn Bk , k  1,.., n Ck , k  2,.., n  A0 y1  y0  B1  A1 x1  x0 B2  B1 y 2  y0  C2  B2 x2  x0 x2  x1 y3  y0 B3  B1  C3  B3 x3  x0 x3  x1 … … y n  y0 Bn  B1  Bn  Cn xn  x0 xn  x1 18 Dk , k  3,.., n  A2 C3  C2  D3 x3  x2 … Cn  C2  Dn xn  x2 … … Берём следующий многочлен B( x)  B( x1 ) C ( x)   A2  A3 ( x  x2 )  ...  An ( x  x2 )...( x  xn1 ). x  x1 B( x2 )  B( x1 ) Тогда С ( x2 )  A2  . x2  x1 И так далее шаг за шагом определяем все Ai , что и даёт нам весь многочлен Ln (x) . Вычисления удобно заносить в таблицу из (n+2)-х столбцов (таблица 1). Пример. Экспериментально получены четыре точки 23,30 24,25 25,25 26,10 xk 299 328 373 415 yk Найдём значение функции y(25) используя все узловые точки, то есть на основе интерполяционного многочлена Лагранжа 3-го порядка. Решение. Заполним разностную таблицу (таблица 2). Таблица 2. Пример расчётной схемы для вычисления многочлена Лагранжа xk Dk , k  3 Ck , k  2,3 Bk , k  1,..,3 yk 23,30 299  A0  299 328  299  30,5263  B1  A1  30,5263 24,25 328 24,25  23,3 373  299  37,9487  30,5263  7,4224  C2  A2  7,4224 25 , 25  23 , 3 25,25 373 25,25  24,25  37,9487  B2 5,8931  7,4224  415  299 41,4286  30,5263   26,1  25,25 26,1  24,25 26,10 415 26,1  23,3  1,7992  D3   41,4286  B3  5,8931  C3  A3  1,7992 Подставляя в формулу (6), получим L3 ( x)  299  30,5263( x  23,30)  7,4224( x  23,30)( x  24,25)   1,7992( x  23,30)( x  24,25)( x  25,25) 19 и y(25)  L3 (25)  299  30,5263 1,70  7,4224 1,70  0,75   1,7992 1,70  0,75  (0,25)  360,9 . Замечание. Формула Лагранжа (3) связывает две переменных x и y, причём любую из них можно считать независимой. Поэтому формулу Лагранжа используют и при задаче обратного интерполирования, то есть нахождения аргумента x по заданному значению функции y с помощью формулы x  Ln ( y) . Пусть теперь интерполяционный многочлен Лагранжа Ln (x) построен для известной функции f (x) . Необходимо выяснить, насколько этот многочлен близок к функции f (x) в точках отрезка [a,b], отличных от узлов интерполяции, то есть оценить погрешность интерполяции | f ( x)  Ln ( x) | . Оценку погрешности можно провести на основании теоремы. Теорема. Пусть функция f (x) (n+1) раз дифференцируема на отрезке [a,b], содержащем узлы интерполяции x0 ,..., xn . Тогда справедливо неравенство: M n 1 (7) | f ( x)  Ln ( x) |  max | Wn 1 ( x) |, (n  1)! [0,6] где M n 1  max | f ( n 1) ( x) |, Wn 1 ( x)  ( x  x0 )( x  x1 )...( x  xn ) . [ 0, 6 ] Для максимальной погрешности интерполяции на всем отрезке [a,b] справедлива оценка M n 1 M n 1 max | f ( x)  Ln ( x) |   max | Wn 1 ( x) |  (b  a) n 1 . [ 0, 6 ] (n  1)! [0,6] (n  1)! Замечание. При n   правая часть оценки стремится к нулю лишь для очень узкого класса функций. Для многих других это не так. Более того, если узлы интерполяции равноотстоят друг от друга, M n 1 то (b  a) n 1   . Если отказаться от условия равноотстоящих (n  1)! углов, то процесс начинает сходиться к породившей его функции. 3.3. Интерполяционный многочлен Ньютона для равноотстоящих узлов Часто интерполирование ведётся для функций, заданных таблично с равноотстоящими узлами. В этом случае шаг h  xi 1  xi , 20 i  0,..., n  1 является величиной постоянной. Для таких таблиц построение интерполяционных формул заметно упрощается. Пусть функция y  f (x) задана таблицей с постоянным шагом x1 y1 x f (x) x2 y2 … … xn yn Разности между значениями функции в соседних узлах интерполяции называется конечными разностями первого порядка yi  yi 1  yi , i  0, 1,..., n  1. Из конечных разностей первого порядка образуются конечные разности второго порядка 2 yi  yi 1  yi i  0, 1,..., n  2 . Таким образом, конечные разности m-го порядка определяются следующим образом m yi  m 1 yi 1  m 1 yi i  0, 1,..., n  m . Конечные разности любого порядка могут быть представлены через значения функции 2 y0  y1  y0  y2  y1  ( y1  y0 )  y2  2 y1  y0 , 2 y1  y2  y1  y3  y2  ( y2  y1 )  y3  2 y2  y1 , …, 3 y0  2 y1  2 y0  y3  2 y2  y1  ( y2  2 y1  y0 )  y3  3 y2  y1 и так далее. Методом математической индукции можно доказать, что m m! m .  yi   (1) k Cmk ym  i  k , где C mk  k ! ( m  k )! k 0 Пример. Выпишем конечные разности различных порядков для функции, заданной таблично (таблица 3). Таблица 3. Разностная схема i 1 2 3 4 5 xi 1 2 3 4 5 yi 0,01 0,02 0,03 0,04 0,05 yi 102 2 yi 102 3 yi 102 4 yi 102 5 yi 102 1 7 19 37 61 6 12 18 24 21 6 6 6 Замечание. Если конечные разности известны, то узлы можно m определить с помощью формулы y m   C mk k y 0 . k 0 Замечание. Если функция f (x) n раз непрерывнодифференцируема на отрезке [a,b] и шаг h достаточно мал (h<<1), то для h-ой производной этой функции имеет место приближенное равенство n f ( x) ( n) , k=0,1,…,n. f ( x)  k h Пусть функция f (x) задана таблично с равноотстоящими узлами. Как и в предыдущем случае, будем искать интерполяционный многочлен в виде: N n ( x)  A0  A1 ( x  x0 )  A2 ( x  x0 )( x  x1 )  ...  An ( x  x0 )...( x  xn1 ) (8) Для определения коэффициентов A0 ,..., An воспользуемся условием совпадения значений исходной функции и многочлена в узлах, т.е. N n ( xi )  yi , i  0, 1,..., n . Тогда N n ( x0 )  A0  y0 , N n ( x1 )  A0  A1 ( x1  x0 )  hA1  A0  y1, N n ( x2 )  A0  A1 ( x2  x0 )  A2 ( x2  x0 )( x2  x1 )   A0  2hA1  A2 2h * h  A0  2hA1  2h 2 A2  y2 , и так далее. Отсюда A0  y0 , A1  y1  A0 y1  y0 y0   , h h h y 2  y0  2 h y1  y0 h   ( A0  2hA1 )  y2  (9) 2h 2 2h 2 y2  y0  2 y1  2 y0 y2  2 y1  y0 y02    2 , …, 2h 2 2h 2 2h k y0 Ak  (методом математической индукции). k k!h Подставляя формулы (9) в (8), получим интерполяционный многочлен Ньютона A2  22 y0 y02 yon N n ( x)  y0  ( x  x0 )  2 ( x  x0 )( x  x1 )  ...  ( x  x0 )...( x  xn1 ), h h h!h n (10) Формулу (10) можно записывать в виде, использующем новую x  x0 переменную t  , то есть x  t  h  x0 и h x  x1 t  h  x0  x1 t  h  h    t  1, h h h x  x2 t  h  x0  x2 t  h  2h    t  2, …, h h h x  xn1  t  n  1. h Таким образом, формулу (10) можно записать в виде: N n ( x)  N n ( x0  th)  (10') y02 y03 y0n  y0  y0  t  t (t  1)  t (t  1)(t  2)  ...  t (t  1)...(t  n  1) 2! 3! n! n или N ( x0  t  h)  y0   Ckt k y0 . k 1 Формула (10') называется первой интерполяционной формулой Ньютона. Эта формула применяется для интерполирования в начале отрезка интерполяции, когда t малό по абсолютной величине (|t|<1). Поэтому формулу (10') называют формулой интерполирования вперёд. С целью повышения точности расчётов и уменьшения числа слагаемых в формуле (10') за начальное приближение x0 можно принимать любое табличное значение xi , при этом формула (10') принимает вид n i N ni ( xi  t  h)  yi   Ckt k yi . k 1 Когда значение аргумента находится ближе к концу отрезка интерполяции применять первую интерполяционную формулу станоx  xn виться невыгодно. В этом случае используется замена    0. h Тогда интерполяционный многочлен Ньютона принимает вид 23 N n ( x)  N n ( xn    h)  (10'') n y0 2 yn2  yn  yn1    (  1)  ...  (  1)...(   n  1) 2! n! n или N n ( xn    h)  yn   Ck k yn  k , и называется вторым интерпоk 1 ляционным многочленом Ньютона или интерполированием назад. Пример. Найдём интерполяционный многочлен Ньютона для функции, заданной таблично, и с его помощью вычислим f (0,596) x f (x) 0,40 0,55 0,70 0,85 1,00 1,15 0,41075 0,57815 0,69675 0,88811 1,02652 1,25386 Решение. Составим таблицу конечных разностей (таблица 4) для интерполирования вперёд, так как x = 0,596 находится ближе к началу отрезка [0,55; 0,70]. Таблица 4. Разностная схема i 1 2 3 4 5 xi 0,40 0,55 0,70 0,85 1,00 1,15 yi yi 105 2 yi 105 3 yi 105 4 yi 105 5 yi 105 0,41075 16740 0,57815 -4880 11860 12156 0,69675 7276 -24727 19136 -12571 51486 -5295 26759 0,88811 13841 14188 8893 1,02652 22734 1,25386 Подставляя в формулу (10'), получаем N 4 ( x)  N 4 (0,55  t  0,15)  0,57815  0,11860  t  0,07276 t (t  1)(t  2)  2! 0,26759 t (t  1)(t  2)(t  3)  0,57815  0,11860t  0,03638t (t  1)  4!  0,02095t (t  1)(t  2)  0,01115t (t  1)(t  2)(t  3) ,  0,96  0,55  f (0,596)  N 4 (0,596)  t    0,57815  0,11860  0,3067  , 15    0,03638  0,3067  (0,3067  1)  0,02095  0,3067  (0,3067  1)   (0,3067  2)  0,01115 * 0,3067(0,3067  1)(0,3067  2)  0,5884 . Отметим, что здесь имеет место нечастый случай, когда разности yk не убывают по абсолютной величине, поэтому формула (10")  24 при i=1 даёт результат, отличный от результата по формуле (10). Действительно, 0,16740 0,04880 N5 ( x)  0,41075  ( x  0,40)  ( x  0,40)( x  0,55)  2 0,15 2  0,15 0,12156  ( x  0,40)( x  0,55)( x  0,70)  3 6  0,15 0,24727  ( x  0,40)( x  0,55)( x  0,70)( x  0,85)  4 24  0,15 0,51486  ( x  0,40)( x  0,55)( x  0,70)( x  0,85)( x  1,0) . 5 120  0,15 f (0,596)  N5 (0,596)  0,41075  1,116  0,196  1,08444  0,00902   6,00296  (0,00094)  20,3514  0,00024  56,50041 (0,00009)   0,6038 . Однако первый способ даёт более точный результат. Оценим погрешность интерполяционного многочлена Ньютона. x  x0 Замена  t позволяет записать оценку h f ( n1) () Rn ( x)  n1 ( x) , где Rn ( x)  f ( x)  N n ( x) , (n  1)! f ( n1) () в виде Rn ( x0  t  h)  (t  h)( x0  t  h  x1 )...( x0  t  h  xn )  n! f ( n1) () h n1 t (t  1)..(t  n)  f ( n1) ()  Cnt 1  h n1 . (n  1)! (n  1)! При достаточно малых h Rn ( x)  Cnt 1y0n1 . Лемма. При  t [0;1] и  n N имеет место неравенство n! 1 , t (t  1)...(t  n)   | Cnt 1 | 4 4(n  1) то есть | Rn ( x) | max{| f ( n1) ()  Cnt 1  h n1 |} , следовательно, | Rn ( x) |  max {| f h n1  M n1 . 4(n  1) почти постоянны и шаг h достаточно мал, ( n 1) x0  xn Если разности y 0n1 () |}  max {| Cnt 1 |}  h n1 0t 1 25 то M n 1h n 1  y0n 1 n 1 y0 .  | Rn ( x) | 4(n  1) (11) 2 y0 При n=1 (линейная интерполяция) | R1 ( x) | . Следователь8 но, если h достаточно мало, а вторые разности почти постоянны и не 1 превышают четырёх единиц младшего разряда, то | R1 ( x) | единиц 2 младшего разряда. Значит, такая линейная интерполяция не оказывает влияние на верные в узком смысле десятичные знаки вычисляемого числа. 3 y0 При n=2 (квадратичная интерполяция) | R2 ( x) | . Однако 12 1 эта оценка сильно завышена. На самом деле | R2 ( x) | единиц 2 младшего разряда, если | 3 y0 | 7 единиц младшего разряда. Таким образом, квадратичная интерполяция допустима, если h достаточно мало, а третьи разности почти постоянны и не превышают 7 единиц младшего разряда. Пример. На основании данных таблицы 4 2 yk почти постоянные, h – мало, 2 yk  1 101  105 единицы младшего разряда (1 105 ), следовательно, линейная интерполяция не обеспечивает пяти верных знаков. 3.4. Аппроксимация функции. Метод наименьших квадратов На практике часто возникает необходимость описать в виде функциональной зависимости связь между величинами, заданными таблично или в виде набора точек с координатами ( xi , yi ) , i=0,…,n. Аппроксимация – это нахождение функции F(x), значения которой на отрезке [ x0 , xn ] как можно меньше отличается от значений аппроксимируемой функции f (x) (рис. 4). Требование точного совпадения, то есть F ( xi )  f ( xi )  yi , приводит к задаче интерполирования, которая уже обсуждалась. 26 Рисунок 4. Выбор аппроксимирующей функции При аппроксимации желательно получить относительно простую функциональную зависимость, например, многочлен, которая позволила бы сгладить экспериментальные погрешности и вычислить значения функции в точках, не содержащихся в исходной таблице. Эта функциональная зависимость F(x) должна с достаточной точностью соответствовать исходной функции f (x) . В качестве критерия точности чаще всего используется метод наименьших квадратов: определяется такая функциональная зависимость F(x), при которой сумма квадратов отклонений минимальна, то есть n S   ( yi  F ( xi )) 2  min . i 0 В 1795 году метод наименьших квадратов был предложен Гауссом и до сих пор является одним из важнейших методов теории ошибок для оценки неизвестных величин по результатам измерений. Погрешность приближения оценивается величиной среднего квадратического отклонения S .  n 1 В качестве функциональной зависимости F(x) чаще всего используется многочлен Pm ( x)  a0  a1 x  a2 x 2  ...  am x m . Тогда условия минимальности суммы квадратов S можно записать, приравнивая нулю частные производные Pm (x) по всем переменным a0 , a1,..., an . Получается система уравнений n P  0   2 ( yi  a0  a1xi  a2 xi2  ...  am xim ) xik  0 k  0,..., m  ak i 0 27 n  (yi  a0  a1xi  a2 xi2  ...  am xim ) xik  0 k  0,..., m  i 0 n a0  i 0 xik n  a1  i 0 xik 1 n  a2  i 0 xik  2 n  ...  am  i 0 n xik  m Если ввести обозначения ck   i 0 xik , n   yi  xik k  0,..., m . i 0 n bk   yi xik , то система i 0 примет вид c0 a0  c1a1  c2 a2  ...  cm am  b0 , c a  c a  c a  ...  c a  b , 1 0 2 1 3 2 n 1 m 1 (12)  ...  cm a0  cm1a1  cm 2 a2  ...  cm m am  bm или в матричной форме записи C  A  B , где С – матрица Грама, симметричная и положительно определённая, основная матрица системы, B – вектор-столбец свободных членов, А – вектор-столбец неизвестных. Рассмотрим частные случаи. Линейная аппроксимация, m = 1. P1( x)  a0  a1x . Элементы основной матрицы системы и вектор-столбца свободных членов n c0   i 0 xi0 n n i 0 n i 0  1  n  1; c1   xi1 n n i 0 i 0 b0   i 0 yi xi0   yi ; b1   n n i 0 i 0   xi ; c2   xi2 ; yi xi1 n   yi xi ; i 0 n    n   n  1  xi    yi     i 0  i 0 C  n n , B   n . 2   xi  xi    yi xi   i 0   i 0  i 0 В силу теоремы Крамера решение системы (12) имеет вид b c b c b c b c a0  0 2 1 21 , a1  1 0 0 21 . c0c2  c1 c0c2  c1 Величина погрешности 1  n 1  n  1 ( yi  (a0  a1 xi )) 2 . i 0 28 Квадратичная аппроксимация, m = 2. P2 ( x)  a0  a1 x  a2 x 2 . Элементы основной матрицы системы и вектор-столбца свободных членов n c0   i 0 xi0 n n i 0 i 0  1  n  1; c1   n c3   i 0 xi3 ; n b1   i 0 n c4   yi xi1 xi4 ; i 0 n xi1 n n i 0 i 0   xi ; c2   xi2 ; n b0   i 0 yi xi0 n   yi ; i 0 n   yi xi ; b2  yi xi2 ; i 0 i 0 n n   n  2  n  1  xi  xi    yi     i 0  i 0 i 0 n n n n     2 3 C    xi  xi  xi , B    yi xi  .  i 0   i 0  i 0 i 0  n 2 n 3 n 4  n  2   xi  xi  xi    yi xi   i 0   i 0  i 0 i 0 В силу теоремы Крамера решение системы (12) имеет вид |C | |C | |C | a0  1 , a1  2 , a2  3 , |C | |C | |C | где C , C1 , C2 , C3 – определители основной и вспомогательных матриц системы. n 1 Величина погрешности  2   ( yi  (a0  a1 xi  a 2 xi2 )) 2 . i 0 n  1 Пример. Построим по методу наименьших квадратов аппроксимирующие многочлены первой и второй степени и оценим погрешность приближения на основе экспериментальных данных 3 5 xi yi 4 6 6 9 7 12 Вычислим коэффициенты c j , b j . Для этого составим таблицу (таб.5) Таблица 5. Расчёты для метода наименьших квадратов xi yi xi 2 xi 3 29 xi 4 yi xi yi xi 2 Сумма 3 5 9 27 81 15 45 4 6 16 64 256 24 96 6 9 36 216 1296 54 324 7 12 49 343 2401 84 588 20 32 110 650 4034 177 1053 Следовательно, коэффициенты c0  4; c1  20; c2  110; c3  650; c4  4034; b0  32; b1  177; b2  1053 . Линейная аппроксимация (m=1) y  a0  a1 x . 32 110  177  20 177  4  32  20 a0    , 5 ; a   1,7  y  0,5  1,7 x . 1 2 2 4 110  20 4 110  20 Квадратичная аппроксимация (m=2) y  a0  a1 x  a2 x 2 . a0  588 2520 120 2  1,63 , a2   7 , a1    0,33  y  7  1,63x  0,33x . 360 360 360 Рисунок 5. Сравнение экспериментальных данных с их аппроксимацией Сравним значения, рассчитанные для линейной и квадратичной аппроксимации, с исходными данными, результаты занесём в таблицу (таб. 6) и изобразим графически (рис. 5). Погрешность приближений составляет 30 1 1,10 ( yi  P1 ( xi )) 2   0,5244 ,  4 4 1 0,10 2  ( yi  P2 ( xi )) 2   0,1581 .  4 4 1  Таблица 6. Сравнение результатов xi 3 4 6 7 yi 5 6 9 12 zi  a0  a1x1 ( yi  zi ) 2 vi  a0  a1x1  a2 x 2 ( yi  vi ) 2 4,6 0,16 5,1 0,01 6,3 0,09 5,8 0,04 9,7 0,49 9,2 0,04 11,4 0,36 11,9 0,01 Сумма= 1,10 Сумма= 0,10 Таким образом, квадратичная функция более точно описывает исходные данные. 3.5. Нелинейная аппроксимация Во многих задачах зависимость между переменными x и y нелинейна по параметрам. Например, y  a  xb – простая модель, но она не является линейной по параметрам a и b. Часто существует преобразование переменных, которое приводит к линейной модели по параметрам. Такое преобразование называется линеаризацией модели. Пример. Линеаризуем математическую модель y  ax b . Решение. Прологарифмируем равенство ln y  ln( ax b )  ln y  ln a  b ln x . Воспользуемся заменой переменных X  ln x, Y  ln y и получим Y  a0  a1 X , где a0  ln a, a1  b . Для ряда часто встречающихся двухпараметрических зависимостей возможные замены, приводящие к линеаризации, представлены в таблице (таб. 7). 31 Вид зависимости Степенная функция y  ax b Показательная функция 1) y  ae bx 2) y  b ae x Гиперболическая функция b 1) y  a  x 2) y  1 a  bx Замена X  ln x Y  ln y X x Y  ln y 1 x Y  ln y X  1 x Yy X  X x 1 Y y Таблица 7 Обратная Примерный график Ограни- замена аппроксимирующей чение для пафункции раметров a  e a0 x>0 ay>0 b  a1 a  e a2 аy>0 b  a1 a  e a2 x  0, аy>0 b  a1 a  a0 x0 b  a1 y0 x a  a0 b a b  a1 32 Логарифмическая функция y  a  b ln x Комбинированная функция 1 y a  b x X  ln x Yy a  a0 x>0 b  a1 X  e x Y 1 y a  a0 y0 b  a1 4. Решение нелинейных уравнений 4.1 Постановка задачи Пусть дана некоторая функция f (x) и требуется найти все или некоторые значения x, для которых (13). f ( x)  0 Значение  , при котором f ()  0 , называется корнем (или решением) уравнения (13). В дальнейшем будем считать что  R (a не С), и f (x) дважды непрерывно-дифференцируема в окрестности корня. Корень  , уравнения (13) называются простым, если первая производная функции f (x) в точке  не равна нулю, т.е. f ()  0 . Если же f ()  0 , то  называется кратным корнем. Рисунок 6. Локализация корней 33 Геометрически корень уравнения (13) – это точка пересечения графика функции y = f(x) с осью абсцисс (рис. 6). Большинство методов отыскания решений уравнения (13) ориентировано на поиск простых корней. 1.2. Основные этапы отыскания решения В процессе приближённого отыскания корней уравнения (13) выделяют два этапа: 1. локализация (или отделение) корней; 2. уточнение корня. Локализация (или отделение) корней – это разбиение области допустимых значений уравнения на интервалы, в каждом из которых содержится не более одного корня. Не существует универсального способа локализации корня. В некоторых случаях отрезок локализации может быть найден из физических соображений. Иногда бывает удобно локализовать корень с помощью графика функции y  f (x) или таблицы её значений. На наличие корня в последнем случае указывает различия знаков функций на концах отрезка (теорема Больцано-Коши). Однако корень чётной кратности таким образом локализовать нельзя (знаки будут одинаковыми). На этапе уточнения корня вычисляют приближенное значение корня с заданной точностью ε > 0. Приближенное значение корня уточняют с помощью различных итерационных методов. Суть этих методов состоит в последовательном вычислении значений x0 , x1 ,..., xn , которые являются приближениями к корню  . 1.3. Метод половинного деления (метод бисекции) Метод деления отрезка пополам является самым простым и надёжным способом решения нелинейного уравнения. Пусть из предварительного анализа известно, что в уравнении (13) корень  [a, b] , причём отрезок этот таков, что f(a)f(b)<0. Разделим отрезок [a,b] на две равные части и получим точку ab и вычислим значения функции в точке x1 . Если f ( x1 )  0 , x1  2 то x1 – корень уравнения (13) и задача решена. Если же f ( x1 )  0 , то рассматриваются два отрезка [a1 , x1 ] и [ x1, b] , причём функция f (x) на концах хотя бы одного из этих отрезков имеет разные знаки. 34 Обозначим такой отрезок [a1 , b1 ] . Очевидно, что корень ba  [a1, b1 ], и длина отрезка b1  a1  . Если b1  a1   , то любую 2 точку отрезка [a1 , b1 ] , (чаще всего середину) можно взять за приближённое значение корня. Если же b1  a1   , то процесс повторяется до тех пор, пока на n-ом шаге не получится либо точное значение корня a b (  xn  n1 n1 ) , либо получится, что bn  an   , и тогда в каче2 b  an стве  можно брать x n , то есть   xn с точностью n . 2 Замечание. Метод половинного деления пригоден для любого уравнения, так как lim xn   , однако сходимость метода невысока. n Она соответствует скорости геометрической прогрессии со знамена1 телем q  . При этом количество итераций можно определить зара2 log 2 (b  a) нее по формуле n   1 , то есть оценка погрешности явля ется априорной. Чтобы увеличить скорость сходимости необходимо использовать малый начальный отрезок [a,b] или, чаще, с помощью метода половинного деления получают малый отрезок, содержащий корень, а затем уточняют корень с помощью других методов. 4.4 Метод касательных Ньютона Метод Ньютона является наиболее эффективным методом решения нелинейных уравнений. Пусть корень  локализирован, то есть  [a, b] , и на отрезке [a,b] других корней нет, пусть f (a)  f (b)  0 , f (x) непрерывна на [a,b] и дважды непрерывно-дифференцируема на интервале (a,b), и пусть первая и вторая производные сохраняют свои знаки на [a,b]. Проведём касательную к графику функции y  f (x) в том из концов отрезка, в котором совпадают знаки функции и её второй производной, то есть f  f   0 . Пусть, например, f (b)  f (b)  0 . Обозначим x0  b – начальное приближение, B0 ( x0 , f ( x0 )) – точка графика y  f (x) , к которой проводится касательная (рис. 7). Уравнения касательной имеет вид 35 y  f ( x0 )  f ' ( x0 )( x  x0 ) . Первое приближение получим, взяв абсциссу точки пересечения f ( x0 ) этой касательной с осью Ox, т.е. x1  x0  . f ' ( x0 ) Аналогично поступим с точкой B1 ( x1, f ( x1 )) , затем с точкой B2 ( x2 , f ( x2 )) и так далее. В результате получим последовательность приближений x0 , x1,...xn ,…, причём f ( xn1 ) . (14) xn  xn1  f ' ( xn1 ) y x2 x1 x0 x Рисунок 7. Реализация метода Ньютона Замечание. Сходимость метода Ньютона зависит от того, насколько близко к корню выбрано начальное приближение. Неудачный выбор начального приближения может дать расходящуюся последовательность. При выполнении условий, наложенных на функцию f (x), последовательность {xn } сходится к  и для n N выполняются неравенства | f ( xn ) | M | xn   | и | xn   | 2 ( xn  xn1 ) 2 , (15) 2m1 m1 где m1  min {| f ' ( x) |} , M 2  max {| f " ( x) |} . a  x b a  x b Неравенства (15) позволяют оценить достигнутую точность nого приближения, если известны m1 и M 2 . Однако для сложных функций сам поиск m1 и M 2 является трудоёмким процессом. По36 этому на практике часто ограничиваются проверкой неравенства | xn  xn1 |  , считая его левую часть достигнутой точностью. В общем же случае, когда некоторые условия на функцию нарушены, используется более надёжный практический критерий окончания итерационного процесса xn  xn1 2   . xn  2 xn1  xn2 Замечание. При выполнении всех условий на функцию f (x) метод касательных имеет квадратичную сходимость, то есть приводит к примерному удвоению числа верных знаков с каждым новым шагом. Кроме того, метод Ньютона можно обобщить на системы нелинейных уравнений. Недостатками метода является его возможная расходимость и необходимость вычислять производную на каждом шаге. 4.5. Метод секущих (метод хорд) Метод хорд является модификацией метода Ньютона, основанный на замене производной её приближением f ( xn )  f ( xn 1 ) . f ' ( xn )  xn  xn 1 Тогда расчётная формула (14) будет иметь вид ( x  x n2 ) f ( x n1 ) x n  x n1  n1 . (16) f ( x n1 )  f ( x n2 ) Это означает, что касательные заменены секущими или хордами (рис. 8). Метод хорд является двух шаговым методом, так как для вычисления приближения xn необходимо вычислить два предыдущих приближения xn1 и x n2 . Для начала итерационного процесса необходимо знать два начальных приближения ab x0 и x1  . 2 В методе хорд используется более практическая формула нахождения x0 a, f (a) f (c)  0; (b  a) f (a)  . x0  b, f (a) f (c)  0; где c  a  f (b)  f (a)  c, f (c)  0;  37 y x3 x2 x1 x0 x Рисунок 8. Реализация метода секущих Очередное приближение xn получается как точка пересечения секущей, соединяющей точки ( xn1 , f ( xn1 )) и ( xn2 , f ( xn2 ) графика y  f (x), с осью Ox. Метод хорд сходится, если  – простой корень уравнения (13) и в некоторой окрестности этого корня функция f (x) дважды непрерывно-дифференцируема, f ( x)  0 . При этом 5 1 | xn   | C | xn1   | p , p   1,618 . 2 Замечание. Метод хорд сходится медленнее, чем метод Ньютона, однако за счёт уменьшения объёма вычислений метод хорд является более предпочтительным. Критерий окончания итераций в данном методе такой же, как и в методе Ньютона | xn  xn1 |  . Пример. Найдём корень уравнения x 3  3x  5  0 тремя методами с точностью ε = 0,001. Построим примерный график функции f ( x)  x 3  3x  5 для того, чтобы локализовать корень (рис. 9). Из графика видно, что  [2;1] , причём f (2)  f (1)  0 . Следовательно, a  2, b  1. 38 Рисунок 9. Локализация корня уравнения x 3  3x  5  0 1. Метод половинного деления. Таблица 8. Расчётная схема метода половинного деления Итерация n 1 2 3 4 5 6 7 8 9 10 Начало отрезка a -2 -1,5 -1,25 -1,188 -1,156 -1,156 -1,156 -1,156 -1,156 -1,154 Конец отрезка b -1 -1 -1 -1,125 -1,125 -1,125 -1,141 -1,148 -1,152 -1,152 Середина отрезка хn -1,5 -1,25 -1,125 -1,1875 -1,15625 -1,14063 -1,14844 -1,15234 -1,1543 -1,15332 Значение функции на концах отрезка и в середине f(a) f(xn) f(b) -9 -2,875 1 -2,875 -0,70313 1 -0,70313 0,201172 1 -0,70313 -0,23706 0,201 -0,23706 -0,01456 0,201 -0,01456 0,094143 0,201 -0,01456 0,040003 0,094 -0,01456 0,012776 0,04 -0,01456 -0,00088 0,013 -0,00088 0,005953 0,013 Начальное приближение Оценка ε 0,5 0,25 0,125 0,0625 0,0313 0,0156 0,0078 0,0039 0,002 0,001 Промежуточный вывод a1=-1,5 b1=-1 a2=-1,25 b2=-1 a3=-1,25 b3=-1,125 a4=-1,1875 b4=-1,125 a5=-1,15625 b5=-1,125 a6=-1,15625 b6=-1,14063 a7=-1,15625 b7=-1,14844 a8=-1,15625 b8=-1,15234 a9=-1,1543 b9=-1,15243 a0  2, b0  1, x0   2 1  1,5 . 2 Остальные вычисления занесём в таблицу (табл. 8). На 10-ом шаге b10  a10  0,000977    0,001    x10  1,15332 , то есть   1,153  0,001. 2. Метод касательных Ньютона. f ' ( x)  3x 2  3; f ( x)  6 x  f (a)  9, f (a)  12; f (b)  1, f (b)  6  Так как f (a)  f (a)  0 , то x0  a  2 – начальное приближение. Остальные вычисления занесём в таблицу (табл. 9). 39 Таблица 9. Расчётная схема метода Ньютона Итерация, n Предыдущий шаг x n 1 f ( xn1 ) f ' ( xn1 ) 1 2 3 4 -2 -1,4 -1,181 -1,155 -9 -1,191 -0,191 -0,002 15 8,88 7,184858 6,99879 xn  xn1  f ( xn1 ) f ' ( xn1 ) -1,4 -1,19108 -1,15453 -1,15417  | xn  xn 1 |  ( xn  xn 1 ) 2 | xn  2 xn 1  xn  2 | 0,6 0,218919 0,026555 0,000354 0,008012 0,000143 0,00000003 На 4-ом шаге   | x4  x3 | 0,000354  0,0001    x4  1,154172 , то есть   1,154  0,001. 3. Метод секущих. Итерация, п Таблица 10. Расчётная схема метода секущих x n2 2 3 4 5 6 -2 -1,5 -1,265 -1,171 -1,155 xn1 f ( x n2 ) f ( xn1 ) xn -1,5 -9 -2,875 -1,26531 -1,265 -2,875 -0,82167 -1,17139 -1,171 -0,82167 -0,12149 -1,15509 -1,155 -0,12149 -0,00645 -1,15418 -1,154 -0,00645 -5,5E-05 -1,15317  | xn  xn1 | 0,234694 0,093917 0,016296 0,000914 0,0000078 (b  a) f (a)  1,1; f (b)  f (a) f (a)  9, f (b)  1, f (c)  0,369; ab Так как f (a)  f (c)  0, то x0  a  2 , x1   1,5 – началь2 ное приближение. Остальные вычисления занесём в таблицу (табл. 10). На 5-ом шаге 1  | x5  x4 | 0,000914  0,001, но   | x5  x4 |  0,000914   0,001, 2 то есть 5-ый шаг не даёт возможности получить верное в узком смысле значение для третьей цифры после запятой. Возникает необходимость в 6-ом шаге. ca 40 1   | x6  x5 |  0,0000078  0,001   0,001    x6  1,15417 , 2 то есть   1,154  0,001. 5. Численное интегрирование функции одной переменной 5.1 Постановка задачи Как известно, далеко не все функции интегрируются в квадратурах, то есть их интегралы не выражаются через элементарные функции. Кроме того, даже если удаётся выразить первообразную F (x) функции f (x) через элементарные функции, она может оказаться сложной для вычислений. Также функция f (x) может задаваться таблично, а значит, точное значение первообразной в это случае получить нельзя. В каждом из перечисленных случаев удобно воспользоваться методами численного интегрирования. Суть численного интегрирования заключается в том, что подынтегральная функция f (x) заменяется другой приближенной функцией так, чтобы интеграл от последней легко вычислялся. 5.2 Метод прямоугольников Формулу прямоугольников можно получить из геометрической интерпретации интеграла. Геометрической интерпретацией определённого интеграла b a f ( x)dx является площадь криволинейной трапе- ции, ограниченной графиком функции y  f (x), осью абсцисс и прямыми x  a, x  b , если f ( x)  0 x [a, b] (рис.10). Рисунок 10. Криволинейная трапеция 41 ba точn ками деления a  x0  x1  ...  xn  b, xi 1  xi  h, i  0, 1,..., n  1. Площадь криволинейной трапеции приближённо равна площади ступенчатой фигуры (рис. 11). Разобьём отрезок [a,b] на n равных частей длиной h  6 y=f(x) 5 4 3 2 1 1 2 3 4 5 6 7 8 Рисунок 11. Ступенчатая фигура Эта фигура состоит из n прямоугольников. Основание i-го прямоугольника образует отрезок [ xi , xi 1 ] длины h и высотой  x  xi 1  f i  . В результате получаем квадратурную формулу средних 2   прямоугольников b b  a n 1  xi  xi 1  (17) I   f ( x)dx  I np   f  . n 2   i 0 a Иногда используют формулы b I   f ( x)dx  I л np b  a n 1    f ( xi ) n i 0 (17') np ba n    f ( xi ) . n i 1 (17") a и b I   f ( x)dx  I п a (17') и (17") называются формулами прямоугольников по левым и правым концам соответственно. Замечание. Точное значение интервала лежит между значениями, полученным формулам (17') и (17"). 42 Для оценки погрешности формулы прямоугольников используется следующая теорема. Теорема. Пусть функция f(x) дважды непрерывнодифференцируема на отрезке [a,b]. Тогда для формулы прямоугольников справедлива оценка M (b  a) 2 | I  I np | 2 h , где M 2  max | f ( x) | . [ a ,b ] 24 5.3. Метод трапеций Формула трапеций выводится так же из геометрических соображений. Рисунок 12. Геометрическая интерпретация метода трапеций Разобьём отрезок [a,b] на n частей точками x0 ,..., xn : a  x0  x1  ...  xn  b. Проведём отрезки прямых x  xi от оси абсцисс до пересечения с графиком функции y  f (x) (рис. 12). Соединим концы отрезков ломаной линией. Тогда площадь криволинейной трапеции приближённо равна площади фигуры, составленной из прямоугольных трапеций. Площадь трапеции, построенной на отрезке [ xi 1, xi ] длины f ( xi )  f ( xi 1 ) ba , равна h  , то получим формулу трапеций h 2 n b I   f ( x)dx  a f ( xn 1 )  f ( xn )   f ( x0 )  f ( x1 ) f ( x1 )  f ( x2 )  I tp  h   ...   2 2 2   43 b  a  f ( x0 )  f ( xn ) n1    f ( xi )  .  n  2  i 1 Теорема (оценка погрешности) Пусть функция f(x) дважды непрерывно-дифференцируема на отрезке [a,b]. Тогда для формулы трапеций справедлива оценка погрешности M (b  a) 2 | I  I тр | 2 h , где M 2  max | f ( x) | . [ a ,b ] 12 Замечание. Сравнение погрешности метода прямоугольников трапеций показывает, что метод прямоугольников в два раза точнее.  5.4. Метод Симпсона (метод парабол) Разобьём отрезок [a,b] на m произвольных частей, а затем каждый у этих m отрезков разобьём на две равные части. В итоге общее число отрезков разбиения n=2m. На каждом из отрезков с чётными концами [ x2k  2 , x2k ] через точки ( x2k 2 , y 2k 2 ) , ( x2k 1 , y 2k 1 ) , ( x2k , y2k ) проведём квадратичную параболу k ( x)  ak x 2  bk x  ck (ось симметрии параллельна оси Оу), то есть воспользуемся квадратичной интерполяцией. В таком случае криволинейная трапеция заменяется набором параболических трапеций (рис. 13). Вычислим площадь одной такой трапеции x2 k x2 k x 2 k 2 x 2 k 2  k ( x)dx   N 2 ( x)dx  x2 k  ( y2k  2  y2k  2  t  x 2 k 2 t (t  1) 2   y2k  2 )dx, 2 x  x2 k  2 в интерполяционном многочлене Ньютона. hk Воспользуемся заменой x  x2 k  2 , x  hk  t  x2k 2  dx  hk dt , tн  0, tв  2, t hk тогда где t  x2 k 2 2 y2k  2 2  k ( x)dx  hk  ( y2k  2  y2k  2  t  2 (t  t ))dt  x 2 k 2 2  y2k  2 2 2 y2k  2  t 2 t 3     hk y2k  2  t  t         2 2  2 3  0  44  2 y2k  2 2  hk   hk   2 y2k  2  2y2k  2     (6 y2k  2  6 y2k 1  2 y2k  2 )  2 3 3  h  k (6 y2k  2  6( y2k 1  y2k  2 )  ( y2k  2 y2k 1  y2k  2 ))  3 h x  x2 k  2 .  k ( y2k  4 y2k 1  y2k  2 ) , где hk  2 k 3 2 m b 1 В результате I   f ( x)dx  I c  hk ( y2k  2  4 y2k 1  y2k ). a 3 k 1 8 7 6 5 4 3 2 1 1 2 3 4 5 6 7 8 Рисунок 13. Геометрическая интерпретация метода парабол ba При разбиении на равные части hk  const  h  получается 2m формула Симпсона b ba m I   f ( x)dx  I c  ( y2k  2  4 y2k 1  y2k ) a 3n k 1 или m m b  ba  y0  yn  4 y2k 1  2 y2k  . I   f ( x)dx  I c  a 6m   k 1 k 1 Теорема (оценка погрешности) Пусть функция f (x) имеет на отрезке [a,b] непрерывную четвертую производную. Тогда для формулы Симпсона справедлива оценка погрешности. M (b  a) 4 | I  I c | n h , где M 4  max | f ( 4) ( x) | . [ a ,b ] 180 Замечание. Формула Симпсона наиболее точная из трех. 45 5.5. Правило Рунге практической оценки погрешности Оценки погрешности формулы прямоугольников, трапеции и Симпсона являются априорными, но для их использования нужно знать M 2 , M 4 , что для сложных функций является утомительным. Оценки погрешности зависят от длины элементарного отрезка h, и при достаточно малом h справедливо приближённое равенство I  I h  Chk . Если уменьшить шаг h в два раза, то h 1 1 I  I 2  k Chk  k ( I  I k ) . 2 2 Вычитая из первого равенства второе, получаем правило Рунге h 2 h h 1 I  Ih k k 2 2 , I  I  k Ch (2  1)  I  I  k 2 2 1 где для формул прямоугольников и трапеций k  2 , для формулы Симпсона k  4 . Это приближённое равенство даёт апостериорную оценку погрешности численного интегрирования. Используя правило Рунге, можно построить процедуру приближённого вычисления интеграла с заданной точностью ε. Нужно, начав вычисления с некоторого шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближённое значение I h . Вычисления заканчиваются, когда выполняется неравенство h | I 2  Ih |  k 2 1 h h h I 2 I 2 и тогда I  I  k . 2 1 1,3 sin x Пример. Вычислим значение интеграла  dx тремя мето0, 3 x дами с шагом h=0,1. sin x Вычислим значения функции f ( x)  на отрезке 0,3; 1,3, x результаты занесём в таблицу 11. Таблица 11. sin x sin x f ( x)  f ( x)  х х x x 0,3 0,985067 0,8 0,896695 0,35 0,979708 0,85 0,883859 46 0,4 0,45 0,5 0,55 0,6 0,65 0,7 0,75 0,973546 0,96659 0,958851 0,95034 0,941071 0,931056 0,920311 0,908852 0,9 0,95 1 1,05 1,1 1,15 1,2 1,25 1,3 0,870363 0,856227 0,841471 0,826117 0,810189 0,793708 0,776699 0,759188 0,741199 1. Метод средних прямоугольников. 1,3 sin x I  dx  I пр  0,3 x 1,3  0,3  f (0,35)  f (0,45)  f (0,55)  f (0,65)  f (0,75)   10  f (0,85)  f (0,95)  f (1,05)  f (1,15)  f (1,25)  0,885565 . Оценим погрешность полученного значения. Имеем cos x sin x   sin x f ( x)    2 2  2 3   f " ( x)  45 . x x   x 45(1,3  0,3) Поэтому | I  I пр |  0,12  0,01875  1,8750 10 2 . 12 2. Метод трапеций. 1,3 sin x I  dx I тр  0,3 x 1,3  0,3  f (0,3)  f (1,3)   f (0,4)  f (0,5)  f (0,6)  f (0,7)   10  2  f (0,8)  f (0,9)  f (1,0)  f (1,1)  f (1,2)  0,885233 . Оценим погрешность полученного значения. Уже известно, что 45(1,3  0,3)  0,12  0,0375  3,75 10 2. f ( x)  45 , поэтому | I  I тр | 12 3. Метод Симпсона. 1,3 sin x I  dx I тр  0,3 x 1,3  0,3  f (0,3)  f (1,3)  4( f (0,35)  f (0,45)  f (0,55)  f (0,75)   6 10  f (0,85)  f (0,95)  f (1,05)  f (1,15)  f (1,25)  2( f (0,4)  f (0,5)   f (0,6)  f (0,7)  f (0,8)  f (0,9)  f (1,0)  f (1,1)  f (1,2)  0,885454 . 47 Для оценки погрешности необходимо оценить производную 4-го порядка. sin x cos x sin x sin x f 4 ( x)   4 2  8 5  24 3  | f 4 ( x) | 2965 x  [0,3;1,3] . x x x x 2965(1,3  0,3) Поэтому | I  I c |  0,14  0,001647  1,647  10  3 . 180 1,3 sin x Пример. Найдём значение интеграла  dx с точностью 0,3 x 3   10 , используя формулу трапеций. Для шага h=0,1 уже получено I h  0,885233 . Уменьшим шаг вдвое h2  0,05 h 1,3  0,3  f (0,3)  f (1,3) I 2  f (0,35)  f (0,4)  f (0,5)  f (0,55)  f (0,6)   20  2  f (0,65)  f (0,7)   f (0,75)  f (0,8)  f (0,85)  f (0,9)  f (0,95)  f (1,0)   f (1,05)  f (1,1)  f (1,15)  f (1,2)  f (1,25))  0,885399 . Согласно правилу Рунге h 1 1  2  для _ трапеции  | I h  I 2 |  | 0,885233  0,885399 | 3 3 3  0,000166  0,166 10  0,001   . Требуемая точность достигнута и h h 2 1,3 sin x h I  I 0,885233  0,885399 I  dx  I 2   0,885399   0,885565 2 0, 3 x 3 2 1 или I  0,886  0,001. 6. Численное решение обыкновенных дифференциальных уравнений 6.1. Постановка задачи Обыкновенное дифференциальное уравнение первого порядка имеет вид y'  f ( x, y( x)) . (18). Решением уравнения (18) является семейство функций y  y(x) , которые при подстановке в уравнение (18) обращают его в верное тождество. График решения дифференциального уравнения называется интегральной кривой (рис. 14). 48 k=tga=f(x,y(x)) 6 5 4 3 2 1 1 2 3 4 5 6 7 8 Рисунок 14. Интегральные кривые Производную y' ( x) в каждой точке (x,y) можно интерпретировать как тангенс угла наклона α касательной к графику решения y  y(x) , проходящего через эту точку. Уравнение (18) определяет целое семейство решений. Чтобы выделить одно решение, используются начальное или граничное условие y( x0 )  y0 , то есть задача Коши y'  f ( x, y( x)), y( x0 )  y0 . (18') В силу теоремы Коши задача Коши имеет единственное решение. Однако даже для простых дифференциальных уравнений аналитическое решение получить не всегда удаётся. Приближённые методы основаны на замене исходного уравнения более простым после обоснованного отбрасывания нескольких содержащихся в нём членов и также не всегда приводят к результату. Численные методы решения обыкновенных дифференциальных уравнений – это алгоритм вычисления приближенных значений функции y  y(x) аргумента x для ряда значений, называемых узлами сетки x0 , x1 ,..., xn , то есть в результате получается таблица значений искомой функции y  y(x). Замечание. Далее к табличным значениям можно применять аппроксимацию или интерполяционный многочлен. Замечание. Обыкновенное дифференциальное уравнение п-го порядка y ( n)  f ( x, y, y' ,..., y ( n 1) ) при помощи замены y ( k ) ( x)  yk ( x), k  0,..., n  2 сводится к системе n обыкновенных дифференциальных уравнений первого порядка. 49 Например, дифференциальное уравнение третьего порядка y  yy" y' x можно свести к системе дифференциальных уравнений первого порядка y0 ' ( x)  y1 ( x)   y1 ' ( x)  y2 ( x)   y ' ( x)  y  y  y  x  2 2 1 Таким образом, достаточно разработать метод решения дифференциальных уравнений первого порядка. 6.2. Метод Эйлера Простейшим методом решения задачи Коши y'  f ( x, y), y0 ( x0 )  y0 является метод Эйлера. В этом случае задача Коши решается на отрезке [ x0 , X ]. Выби- X  x0 и строится сетка с системой узлов xi  x0  ih , n i  0,..., n. В методе Эйлера вычисляются приближенные значения искомой функции y  y(x) в узлах сетки yi  y( xi ) с помощью замены производной y' ( x) конечными разностями на отрезках [ xi 1 , xi ], i  1,..., n, то есть yi  yi 1  f ( xi , yi ), i  1,..., n  h yi  yi 1  f ( x0 , yi 1 )h, i  1,..., n, y0  y( x0 ) . (19) Формулы (19) являются расчётными формулами метода Эйлера. Геометрически один шаг метод Эйлера означает замену на отрезке [ xn 1, xi ] интегральной кривой y  y(x) на её касательную y  y' ( xi )( x  xi )  yi . После выполнения n шагов неизвестная интегральная кривая заменяется ломаной Эйлера. Теорема (оценка погрешности метода). Пусть функция f ( x, y) удовлетворяет на отрезке [ x0 , X ] условиям рается шаг h  f f f f  K,   f  L. y x x y Тогда для метода Эйлера справедлива оценка погрешности 50 l 2 L KL l 2 h KL max y ( xi )  yi  e  e , 0i  n 2n 2 где l  X  x0 , то есть метод Эйлера имеет первый порядок точности. Замечание. Оценка точности метода Эйлера бывает затруднена, так как требует вычисления частных производных функции f(x,y). На практике для процедуры приближённого вычисления решения задачи Коши с точностью ε используют правило Рунге: начинают вычисления с некоторого значения шага h, затем уменьшают шаг в h h два раза и вычисляют приближённые значения yih , yi 2 , i  0,..., n . 2 Если | yih  yih / 2 |   , i  0,..., n , то yih / 2 , i  0,..., n будут приближёнными значениями исходной функции y  y(x) в узлах сетки. Замечание. На данный момент существуют несколько модификаций метода Эйлера, которые обеспечивают уже второй порядок точности. 6.3. Метод Рунге-Кутта Рассмотрим задачу Коши (18'), выберем отрезок [ x0 , X ] и шаг X  x0 , построим сетку с системой узлов xi  x0  ih, i  0,..., n, h n yi  y( xi ), которые получаются на основе расчётных формул метода Рунге-Кутта четвёртого порядка точности h (20) yi  yi 1  (k1  2k2  2k3  k4 ) , 6 h h  где k1  f ( xi 1, yi 1 ), k2  f  xi 1  , yi 1  k1 , 2 2  h h  k3  f  xi 1  , yi 1  k2  , k4  f ( xi  h, yi 1  k3h), i  1,..., n . 2 2  Если функция f ( x, y) непрерывна и ограничена вместе со своими производными до 4-го порядка включительно, то равенство (20) обеспечивает хорошие результаты. Априорная оценка погрешности полученных значений yi очень сложна и на практике не используется. 51 Грубую оценку погрешности даёт правило Рунге (Рунге| yih / 2  yih | Румберга). Если   для всех i=0,..,n, то yih / 2 , i  0,..., n бу15 дут приближенными значениями искомой функции y  y(x) в углах сетки. Метод Рунге-Кутта является одним из наиболее распространённых численных методов высокой точности. Метод Эйлера можно рассматривать как простейший вариант метода Рунге-Кутта. И метод Эйлера, и метод Рунге-Кутта являются одношаговыми, то есть позволяют начать счёт сразу с i=0 по известным начальным условиям. Кроме того, эта особенность допускает изменение шага в любой точке в процессе счета, что позволяет строить числовые алгоритмы с автоматическим выбором шага. Для решения дифференциальных уравнений разработаны более точные и быстрые методы, называемые методами прогноза и коррекции, однако они не являются «самостартующими», то есть они – многошаговые, а значит, для получения исходных данных используется какой-либо одношаговый, чаще всего метод Рунге-Кутта. Коррекция в таких методах необходима для определения следующего шага h к узлу сетки. y 12 Пример. Решим задачу Коши y '   3 , y (1)  8 аналитичеx x ски, методом Эйлера и методом Рунге-Кутта с шагом h=0,1 на отрезке [1,2]. y 12 Здесь f ( x, y )   3 , x0  1, y0  8 , а точное решение задачи x x 4 Коши y  2  4 x . x Метод Эйлера yi  yi 1  f ( x0 , yi 1 )h . h Метод Рунге-Кутта yi  yi 1  (k1i  2k2i  2k3i  k4i ) . 6 Расчёты и узловые точки занесены в таблицу 12. Графически не трудно убедится (рис. 15), что кривая, полученная методом Рунге-Кутта ближе к интегральной кривой по сравнению с кривой, построенной методом Эйлера. Пример. Решим задачу Коши y  6 y  13 y  6 cos x  e 2 x , 52 y(0)  2, y(0)  1 аналитически, методом Эйлера и методом Рунге-Кутта с шагом h=0,1 на отрезке [0,1]. Здесь ( x, y, z )  z, ( x, y, z )  6 cos x  e2 x  6 z  13 y , x0  0 , y0  2 , z0  1, а точное решение задачи Коши e3x y  2e x  2e3 x sin x  31sin 2 x  4e3 x cos x  22 cos 2 x . 10 Метод Эйлера yi  yi 1  f ( x0 , yi 1 )h . h Метод Рунге-Кутта yi  yi 1  (k1i  2k2i  2k3i  k4i ) . 6 Расчёты, выполненные в ППП Excel, полученные узловые точки и графики изображены на рисунке 16.   53 Таблица 12. Расчётная схема методов Эйлера и Рунге-Кутта для уравнения 1-го порядка i xi 1 2 3 4 5 6 7 8 9 10 2 2,5 3 3,5 4 4,5 5 5,5 6 6,5 7 точное y 4 7,8125 13,5 21,4375 32 45,5625 62,5 83,1875 108 137,3125 171,5 Эйлер yi 4 7 11,525 17,94583 26,63452 37,96384 52,30704 70,03775 91,52982 117,1573 147,2944 Рунге yi 4 7,812294 13,49961 21,43694 31,99928 45,56162 62,49897 83,18633 107,9987 137,311 171,4984 k1 xi+h/2 yi+k*h/2 k2 xi+h/2 yi+k*h/2 k3 xi+h yi+k*h k4 6 9,374918 13,49987 18,37484 23,99982 30,37481 37,49979 45,37479 53,99978 63,37478 2,25 2,75 3,25 3,75 4,25 4,75 5,25 5,75 6,25 6,75 5,5 10,15602 16,87458 26,03065 37,99923 53,15532 71,87392 94,53003 121,4986 153,1547 7,506944 11,2556 15,75468 21,00401 27,0035 33,75309 41,25277 49,5025 58,50228 68,25209 2,25 2,75 3,25 3,75 4,25 4,75 5,25 5,75 6,25 6,75 5,876736 10,62619 17,43828 26,68794 38,75015 53,9999 72,81217 95,56195 122,6243 154,3741 7,674383 11,42657 15,92812 21,17928 27,18018 33,9309 41,43148 49,68197 58,68238 68,43273 2,5 3 3,5 4 4,5 5 5,5 6 6,5 7 7,837191 13,52558 21,46367 32,02658 45,58937 62,52707 83,21472 108,0273 137,3399 171,5274 9,384877 13,50853 18,38248 24,00665 30,38097 37,50541 45,37995 54,00455 63,37921 73,50392 54 Рисунок 15. Графическое представление результатов численного решения дифференциального уравнения Рисунок 16. Реализация расчётных схем методов Эйлера и Рунге-Кутта в ППП Excel 55 ЛИТЕРАТУРА 1. Сборник задач по математике для втузов. Ч. 2. Специальные разделы математического анализа : учеб. пособие для инж.-техн. спец. вузов / под ред. А. В. Ефимов ; под ред. Б. П. Демидович . - М. : Наука , 2007. - 368 с., ил. 2. Сборник задач по математике для втузов. Ч. 1. Линейная алгебра : учеб. пособие для инж.-техн. спец. вузов / под ред. А. В. Ефимов ; под ред. Б. П. Демидович . - М. : Наука , 2007. - 366 с., ил. 3. Бахвалов, Н.С. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.Н. Кобельков. - М. : БИНОМ. Лаборатория знаний, 2003 4. Зайцева М. И. Элементы вычислительной математики : Учебное пособие / М. И. Зайцева, В. Н. Макаров, А. А. Сидоренко - М ., 2003.- 147 с. Учебное издание Семенова Галина Александровна Савостикова Татьяна Владимировна Математика. Численные Учебное пособие Редактор Г.В. Каркушина Технический редактор Федеральное государственное образовательное учреждение высшего профессионального образования «Государственный университет – учебно-научно-производственный комплекс» Лицензия ИД № 00670 от 05.01.2000г Подписано к печати Формат 60x84 1/16 Усл.печ. л. Тираж экз. Заказ № Отпечатано с готового оригинал - макета 57