close

Вход

Забыли?

вход по аккаунту

?

Grishanovasokmet

код для вставкиСкачать
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
Государственное образовательное учреждение
высшего профессионального образования
САНКТ$ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
ВВЕДЕНИЕ
В ЧИСЛЕННЫЕ МЕТОДЫ
Методические указания
к выполнению лабораторных работ № 1–3
Санкт$Петербург
2005
Составители: Л. И. Гришанова, М. В. Соколовская
Рецензент канд. техн. наук, доцент В. П. Попов
Под редакцией д$ра техн. наук, проф. Р. И. Сольницева
Методические указания предназначены для первоначального изу$
чения численных методов и алгоритмов решения нелинейных алгеб$
раических уравнений, вычисления определенных интегралов и мето$
дов численного интегрирования обыкновенных дифференциальных
уравнений. Методические указания используются в учебном процессе
при чтении курсов «Модели и методы анализа проектных решений»,
«Моделирование в проектировании и производстве» для студентов
очной и заочной форм обучения специальности 220300.
Подготовлены кафедрой компьютерных систем автоматизации и
рекомендованы к изданию редакционно$издательским советом Санкт$
Петербургского государственного университета аэрокосмического при$
боростроения.
Редактор Г. Д. Бакастова
Компьютерная верстка А. Н. Колешко
Сдано в набор 09.11.04. Подписано к печати 27.04.05. Формат 60´84 1/16.
Бумага офсетная. Печать офсетная. Усл. печ. л. 3,31. Уч. $изд. л. 3,46. Тираж 100 экз.
Заказ №
Редакционно$издательский отдел
Отдел электронных публикаций и библиографии библиотеки
Отдел оперативной полиграфии
СПбГУАП
190000, Санкт$Петербург, ул. Б. Морская, 67
©
2
ГОУ ВПО «СПбГУАП», 2005
ПРЕДИСЛОВИЕ
Данные методические указания предназначены для первоначаль$
ного изучения численных методов и алгоритмов решения нелиней$
ных уравнений, вычисления определенных интегралов, а также ме$
тодов численного интегрирования обыкновенных дифференциальных
уравнений. Первая часть методических указаний содержит три лабо$
раторных работы:
– № 1 – «Численное решение нелинейных уравнений»,
– № 2 – «Вычисление определенных интегралов»,
– № 3 – «Численное интегрирование обыкновенных дифференци$
альных уравнений».
В методических указаниях для каждой работы имеется описание
применяемых методов на простых примерах, а также варианты зада$
ний для их выполнения. В первой части методических указаний опи$
саны три первые лабораторные работы.
3
Лабораторная работа № 1
ЧИСЛЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНЫХ УРАВНЕНИЙ
1.1. Общие методические указания к численным методам
решения нелинейных уравнений
Уравнением называется равенство
f(x) = 0,
(1.1)
справедливое при некоторых значениях x = x*, называемых корнями
этого уравнения или нулями функции f(x). Решение уравнения зак$
лючается в определении его корней. Среди корней x* могут быть и
комплексные, однако в данной работе вычисляются только действи$
тельные корни.
Вычисление каждого из действительных корней складывается из
двух этапов:
1) отделение корня, т. е. нахождение возможно малого интервала
[a, b], в пределах которого находится один и только один корень x*
уравнения;
2) уточнение значения корня, т. е. вычисление с заданной степе$
нью точности.
При использовании рассматриваемых ниже методов решения урав$
нения (1.1) к функции f(x) на интервале [a, b] предъявляются следу$
ющие требования:
а) функция f(x) непрерывна и дважды дифференцируема (т. е. су$
ществует первая и вторая производные);
б) первая производная f'(x) непрерывна, сохраняет знак и не обра$
щается в нуль;
в) вторая производная f"(x) непрерывна и сохраняет знак.
Отделение корней может производиться аналитическим или гра$
фическим способами. Аналитический способ основывается на теоре$
ме Коши, утверждающей, что для непрерывной функции f(x) (первое
требование «а»), принимающей на концах интервала [a, b] разные
знаки, т. е. f(a)×f(b)< 0, уравнение (1.1) имеет внутри этого интерва$
ла хотя бы один корень (рис. 1.1, а). Если к этому добавить второе
требование «б», означающее монотонность функции f(x), то этот ко$
4
рень оказывается единственным (рис. 1.1, б). В противном случае
следует интервал [ak, ak+1] разделить на меньшие интервалы, повто$
ряя для каждого из них указанные действия. В этих условиях отде$
ление корня сводится к вычислению значений функции f(x) для пос$
ледовательности точек a1, a2, …, an и сопоставлению знаков f(ak),
f(ak+1) в соседних точках ak и ak+1.
При использовании графического способа уравнение (1.1) следует
представить в виде
f1(x) = f2(x)
(1.2)
и построить графики функций y = f1(x) и y = f2(x). Абсцисса точки
пересечения этих графиков дает приближенное значение x0 корня x*
уравнения
f(x) = f1(x) – f2(x) = 0.
Представление уравнения (1.1) в форме (1.2) не является, есте$
ственно, однозначным и его следует подбирать так, чтобы построе$
ние графиков было возможно простым.
Из того же чертежа следует определить и тот интервал [a, b], в
пределах которого данный корень является единственным (если это
необходимо для выбранного метода последующего уточнения значе$
ния корня x0).
При необходимости отделения (и вычисления) отрицательного
корня x* достаточно в уравнении (1.1) сделать замену Z = – x. Тогда
каждому отрицательному корню x* исходного уравнения (1.1) соот$
ветствует положительный корень Z* уравнения f(Z) = 0; после на$
хождения этого корня остается положить x* = –Z*.
Пример 1.1
Отделить положительные корни уравнения x3 –x –1 = 0.
a)
1
б)
1
Рис. 1.1. Корни уравнения x3 – x –1 = 0.
5
1
Р е ш е н и е . Выберем y = f1(x) = x3, y =
= f2(x) = x+1.
Построение графиков (рис. 1.2) показыва$
ет, что исходное уравнение имеет единствен$
ный положительный корень x*, заключенный
в интервале [1.0;1.5], приближенное значение
которого можно принять равным x0 = 1.25.
Выделенный корень является в данном случае
единственным положительным корнем урав$
Рис. 1.2. Отделение
нения, так как других точек пересечения гра$
корней уравнения
фиков нет.
После того как корень x* отделен (определен интервал, в котором
он является единственным), проводится уточнение значения корня.
Для решения этой задачи разработано большое число различных ме$
тодов, целесообразность применения каждого из которых определя$
ется особенностями функции f(x).
При выполнении данной работы используются наиболее простые
из них, опирающиеся на свойства функции f(x), определяемые ука$
занными ранее требованиями «а», «б», «в».
1.2. Метод половинного деления
Пусть в уравнении (1.1) функция f(x) является непрерывной (пер$
вое требование «а») на интервале [a, b], в котором расположен один
искомый корень x*. Для нахождения этого корня разделим отрезок
[a, b] пополам точкой
a1b
.
2
Если теперь f(x1) = 0, то x1 и является корнем уравнения. В про$
тивном случае выбираем тот из отрезков [a, x1 ] или x[1, b], на концах
которого функция f(x) имеет разные знаки. В пределах этого отрезка
согласно предыдущим рассуждениям лежит искомый корень. Таким
образом, оказывается определенным интервал [a1, b1] (где a1 = a, b1 = x1
или a1 = x1, b1 = b), меньший первоначального [a, b] и содержащий x*.
Повторяя подобные построения, получаем последовательность умень$
шающихся интервалов [an, bn] таких, что
23
bn 1 an 2
1
(b 1 a)
(1.3)
2n
*
и в каждом из которых заключен корень x . Точность вычисления
*
корня x определяется размерами интервала [an, bn] после n$го деле$
6
ния исходного интервала [a, b], так как ошибка определения корня
x* не превышает величины bn– an.
Следовательно, если e есть заданная точность вычисления, то дол$
жно выполняться условие
1
(b 1 a) 2 3.
(1.4)
2n
Отсюда можно определить и необходимое число шагов половин$
ного деления интервала [a, b], если задано e:
b1a
2 4 3.32lg b 1 a .
lg2
2
lg
n3
Пример 1.2
Методом половинного деле$
ния найти корень уравнения
f(x) = x4+2x3–x–1 = 0, лежащий
в интервале [0, 1], после шести
шагов половинного деления.
Р е ш е н и е . Полагая a = 0,
b = 1, последовательно находим
(рис. 1.3)
Рис. 1.3. Метод половинного
деления
f(0) = – 1.0; f(1) = 1.0;
0 11
3 0.5;
2
f(0.5) = 0.06+0.25 – 0.5 – 1 = – 1.19 < 0.
21 3
Отсюда
f(0.5)×f(0) > 0, f(0.5)×f(1) < 0.
Следовательно, необходимо взять
22 3
0.5 1 1.0
3 0.75;
2
f(0.75) = 0.32+0.84 – 0.75 –1 = 0.59< 0,
f(0.75)×f(1) < 0; a2 = x2 = 0.75; b2 = b1 = 1.0
23 3
0.75 1 1.0
3 0.875;
2
7
f(0.875) = 0.05> 0,
f(.0875)×f(0.75) < 0; a3 = a2 = 0.75; b3 = x3 = 0.875
24 3
0.75 1 0.875
3 0.8125;
2
f(0.8125) = – 0.304 < 0
f(0.8125) ×f(0.75)< 0; a4 = x4 = 0.8125; b4 = b3 = 0.875
25 3
0.8125 1 0.875
3 0.8438;
2
f(0.8438) = – 0.135 > 0;
f(0.8438) ×f(0.8125) < 0; a5 = x5 = 0.8438; b5 = b4 = 0.875
26 3
0.8438 1 0.875
3 0.8594;
2
f(0.8594) = – 0.43 < 0; f(0.8594) ×f(0.875)< 0.
Следовательно, a6 = x6 = 0.8594; b6 = b5 = 0.875.
После шестого деления исходного интервала [0, 1] установлено,
что искомый корень x* удовлетворяет условию 0.8594 < x* < 0.8750
и максимальная погрешность e его вычисления составляет e <
< 0.8750– 0.8594 = 0.0156. Для уменьшения этой погрешности мож$
но принять
1
x* 1 (0.8594 2 0.875) 3 0.867.
2
1.3. Метод последовательных приближений
Для применения этого метода следует заменить исходное уравне$
ние (1.1) равносильным (имеющим те же корни) уравнением вида x =
j(x). (1.5). Пусть теперь каким$либо образом выбрано приближен$
ное значение x0 искомого корня x*. Построим числовую последова$
тельность x1, x2, …, пользуясь выражениями
xn+1 = j(xn) (n = 0, 1, 2, …).
Если функция j(x) непрерывна и существует предел, то
3
4
3
4
lim xn 11 5 lim 1 6(xn ) 2 5 6 7 lim xn 8 5 6 7 lim xn 11 8 .
n 23
9n23
9n23
n 23
8
(1.6)
Таким образом, предел xn+1 является корнем x* уравнения (1.5)
lim xn 11 1 x*, а следовательно, и корнем исходного уравнения (1.1),
n 23
который может быть вычислен по формуле (1.6) с любой степенью
точности. Геометрически такой метод можно пояснить следующим
образом (рис. 1.4, а). Построим графики функций y = x и y = j(x),
точка M пересечения которых определяет искомый корень x*. Начи$
ная с точки A0[x0, j(x0)], после первого вычисления x1 = j(x0) опре$
деляем точку B1[x1, j(x0)] и значение x1 первого приближения. Сле$
дующее вычисление x2 = j(x1) определяет точку B2[x2, j(x1)]. Повто$
рение таких вычислений дает последовательность точек B2, B3,...,
абсциссы x2, x3,... которых представляют собой последовательные
приближения к корню x*.
Приведенное построение отвечает случаю, когда j'(x) > 0 (после$
довательность x0, x1,... сходится к x* монотонно).
Если j'(x) < 0, то процесс приближений не имеет монотонного ха$
рактера (рис. 1.4, б).
a)
б)
в)
Рис. 1.4. Метод последовательных приближений
9
Если для интервала [a, b] выполняется условие
êj'(x) ê£ m < 1,
(1.7)
где m – наибольшее значение êj'(x) ê на этом интервале, то погреш$
ность e n$го приближения может быть оценена неравенством
m*
x11x0 ,
(1.8)
13m
*
и последовательные приближения сходятся к x при любом выборе
x0 на интервале [a, b].
В том случае, когда функция j(x) выбрана в виде
1 2 xn 3 x * 4
j(x) = x – f(x),
(1.9)
*
имеет место соотношение çxn+1 – x ç£ çxn+1 – xnç, из которого следует,
что при заданной точности e приближения вычисления прекращают$
ся при выполнении условия
çxn+1 – xnç £ e.
(1.10)
При выполнении условия (1.7) из соотношения (1.8) следует, что
при увеличении числа n проводимых вычислений величина çxn – x*ç
неограниченно возрастает и последовательные приближения не схо$
дятся к корню x* (рис. 1.4, в).
Пример 1.3
Найти корень уравнения sin x – 2x+0.5 = 0, расположенный в
интервале [0, p/2], ограничиваясь точностью вычислений e £ 0.0001.
Р е ш е н и е . Заменим исходное уравнение равносильным, т. е.
представим его в виде (1.5).
Тогда
x = 0.25+0.5sin x = j(x),
j'(x) = 0.5cos x и | j'(x)| £ 0.5 < 1
при любых x, в частности при 0 £ x £ p/2, и условия сходимости со$
блюдены.
Формула (1.6) метода последовательных приближений принима$
ет вид
xn+1 = 0.25+0.5sin xn (n = 0, 1, 2,...).
Начиная с точки x0 = 0.5 рассматриваемого интервала, последо$
вательно находим следующие значения:
x1 = 0.25+0.5×sin 0.5 = 0.4897,
10
x2 = 0.25+0.5×sin 0.4897 = 0.4852,
x3 = 0.25+0.5×sin 0.4852 = 0.4832,
x4 = 0.25+0.5×sin 0.4832 = 0.4823,
x5 = 0.25+0.5×sin 0.4823 = 0.4819,
x6 = 0.25+0.5×sin 0.4819 = 0.48175,
x7 = 0.25+0.5×sin 0.48175 = 0.48165,
x8 = 0.25+0.5×sin 0.48165 = 0.4816,
x8 » x9 » x10 »ј
Последующие приближения удовлетворяют условию |xn–1 – xn| £ e
= 0.0001 и x* » 0.4816. Его можно принять в качестве значения иско$
мого корня.
1.4. Метод последовательных приближений с переходом
к обратным функциям
Как указывалось при рассмотрении метода последовательных при$
ближений, использование формул (1.6) вида xn+1 = j(xn) (n = 0, 1, 2, …)
требует выбора j(x) с соблюдением условия (1.7), обеспечивающего
сходимость последовательных приближений x0, x1,... к искомому
корню x*. Ранее обращалось внимание на необходимость надлежа$
щего выбора j(x). На этом пути возможным является следующий об$
щий прием.
Если для уравнения (1.5) вида x = j(x) в окрестности искомого
корня (на интервале [a, b]) имеет место êj'(x)ê³ M > 1 и не обеспечи$
вается условие сходимости, то это уравнение следует заменить на
равносильное, получаемое переходом в (1.5) к обратным функци$
ям x = Y(x), где Y(x) – функция, обратная j(x) (функция x обратная
себе самой). Тогда в силу
1 '(x) 2
1
1 '1 (1(x))
имеем
1 '(x) 2
1
1
3
11 ,
1 '(x) M
и процесс последовательных приближений сходится.
Пример 1.4
Требуется найти методом последовательных приближений поло$
жительные корни уравнения 5x – 8(ln x+1) = 0.
11
Представляя это уравнение в
виде
5
x 1 1 2 ln x
8
и строя графики функций y =
= 5/8 x – 1 и y = ln x (рис. 1. 5),
приближенно находим x *1 »
» 0.45, x *2 » 3.7.
Для уточнения правого
корня x*2 » 3.7 методом после$
довательных приближений
Рис. 1.5 Отделение
находим x = 1.6(1+ln x) = j(x),
корней уравнения
и так как в окрестности этого
корня j'(x) = 1.6/x £ 1, то про$
цесс сходится (если x0 выбрано достаточно близко к x *2).
Однако в окрестности корня x*1» 0.45 j'(x)» 3.5 > 1, и процесс
расходится. Поэтому следует перейти к равносильному уравнению
относительно обратных функций x = е0.625x–1 = Y(x).
Как видно, в окрестности корня Y'(x) = 0.625е0.625x–1 » 0.3 < 1, и
последовательные приближения сходятся.
1.5. Метод касательных (метод Ньютона)
Пусть определен интервал [a, b], в пределах которого расположен
единственный корень x* уравнения (1) f(x) = 0, и функция f(x) удов$
летворяет требованиям «а», «б», «в», (см. п. 1.1). Пусть в качестве$
начального приближения выбрано значение x0 и отвечающая ему
точка A0 [x0, f(x0)] (рис. 1.6).
Найдем следующее приближение как точку пересечения с осью,
касательной к кривой y = f(x), проведенной в точке A0. Определив
таким образом новую точку A1 [x1, f(x1)], найдем следующее прибли$
жение x2 как абсциссу пересечения касательной в точке A1.
Повторяя подобные операции, построим последовательность x1,
x2,..., сходящуюся (при выполнении определенных условий) к иско$
мому корню x*.
Для некоторого приближения xn уравнение касательной, проведен$
ной в точке An [xn, f(xn)], имеет вид y – f (xn) = f'(xn)(x – xn), и значение
xn+1, отвечающее следующему приближению, находится из условия
y(xn+1) = 0, отвечающего пересечению касательной с осью абсцисс. Отсюда
xn 11 1 xn 2
12
f (xn )
(n 1 0,1,2,....),
f '(xn )
(1.11)
что и определяет последовательность расчетных формул метода.
Можно показать, что если f(x) удовлетворяет упомянутым требова$
ниям «а», «б», «в» и начальное приближение выбрано так, что
f(x0)×f"(x0)>0,
(1.12)
где f"(x0) = f"(x)½x = x0, то последовательные приближения x0, x1,...
*
сходятся к единственному на интервале [a, b] корню x , который
может быть вычислен с любой необходимой точностью.
Точность вычислений может быть оценена с помощью соотношения
|xn+1 – x* | £ | xn+1 – xn|,
(1.13)
из которого следует, что вычисления могут быть прекращены с получе$
*
нием заданной точности | xn+1– x |£ e, когда выполняется неравенство
| xn+1 – xn|£ e,
(1.14)
аналогичное условию (1.10) метода последовательных приближений.
При несоблюдении условия (1.12) выбора начальной точки x0 ме$
тод может дать неудовлетворительный результат. Например, при
выборе x0 = a (рис. 1.6) последующее
приближение x'1 может оказаться вне
интервала [a, b] отделения данного кор$
ня x*. При этом в соответствии с (1.12)
в качестве начальной следует выбрать
ту из граничных точек, в которой знаки
функции и ее второй производной совпа$
дают.
Пример 1.5
Вычислить наименьший положи$
Рис.1.6. Метод касательных
тельный корень уравнения еx–3x = 0 с
четырьмя десятичными знаками (e =
= 0.0001).
Р е ш е н и е . Для отделения возмож$
ных положительных корней перепишем
уравнение в виде еx = 3x и построим гра$
фики функций y = еx, y = 3x (рис. 1.7),
показывающие, что имеются два поло$
жительных корня, наименьший из ко$
торых x*1 заключен в интервале [0, 1].
Проверяя выполнение требований
«а», «б», «в», предъявляемых к функ$
Рис. 1.7. Отделение
ции f(x) = еx – 3x, находим
корней уравнения
13
f'(x) = еx – 3, f"(x) = еx,
f'(x)<0, f"(x)>0 при 0 £ x £ 1
производные непрерывны и сохраняют знаки на рассматриваемом
интервале. При этом f(x) так же непрерывна, но f(0)× f(1)<0.
Выбирая начальную точку x0, находим
f(0)× f"(0) = (еx– 3x) еxЅx = 0 = 1>0,
f(1)× f"(1) = (еx– 3x) еxЅx = 1 = (е – 3)е < 0.
Следовательно, необходимо принять x0 = 0. Последующие при$
ближения вычисляются по формулам (1.11), которые в данном слу$
чае принимают вид
xn 11 2 xn 1
exn 1 3xn
(n = 0, 1, 2, …).
exn 1 3
Результаты следующих отсюда вычислений сведены в табл. 1.1.
f (xn )
2 xn 11 3 xn в соответствии с (1.13) и
f '(xn )
(1.14) определяет точность приближения xn+1 к искомому корню x*.
Таблица 1.1
При этом величина 1n 2
n
xn
ex
3x
f(xn)
f'(xn)
en
0
0.00000
1.0000
0.0000
1.0000
–2.0
–0.500000
1
0.50000
1.6500
1.5000
0.1500
–1.35
–0.110000
2
0.61000
1.8404
1.8300
0.0104
–1.16
–0.009000
3
0.61900
1.8570
1.8570
0.0001
–1.14
–0.000088
4
0.61909
1.8573
1.8573
0.0000
–
–
В качестве искомого значения корня принимаем x » 0.619.
1.6. Упрощенный метод касательных
Как видно из формул (1.11), их использование требует вычисле$
ния значения производной f'(xn) в каждой из точек xn. Метод можно
упростить, сохраняя для всех последовательных значений x0, x1,...
одно и то же значение производной f'(xn), равное ее значению в на$
чальной точке x0, т. е. полагая
f'(xn)» f'(x0).
14
(1.15)
Тогда формулы (1.11) метода
касательных принимают вид
xn 11 1 xn 2
f (xn )
f '(x0 )
(n 1 0,1,2,...).
(1.16)
Графическое представление та$
кого метода приведено на рис. 1.8,
на котором угол наклона всех ка$
Рис. 1.8. Упрощенный метод
сательных одинаков и определяет$
касательных
ся значением f'(x0). Все остальные
особенности метода касательных сохраняются; прекращение вычис$
лений определяется условием (1.14).
1.7. Метод хорд
Пусть определен интервал [a, b], в котором лежит один корень x*
уравнения (1.1) f(x) = 0.
Учитывая, что f(a)×f(b)<0, определяем первое приближение как
точку пересечения с осью абсцисс хорды A0B0, соединяющей точки
A0[a, f(a)] и B0[b, f(b)] (рис. 1.9, а).
Для нахождения последующего приближения вычислим значе$
ние f(x1) и сопоставим со значениями f(a) и f(b). Выберем тот из ин$
тервалов [a, x1] или [x1, b], на концах которого функция f(x) имеет
разные знаки (именно внутри этого интервала лежит искомый ко$
рень x*). Применим предыдущий прием к этому интервалу, получая
последующее приближение – точку x2.
Заметим, что для случая, показанного на рис. 1.9, а, производные
f'(x) и f"(x) сохраняют положительный знак (f'(x)>0, f"(x)>0;
f'(x)×f"(x)>0) и все приближения x1, x2,... образуют возрастающую
последовательность, ограниченную значением x = x*. Следователь$
но, lim xn 1 x * и при этом в любом из приближений соответствую$
n
12
щая хорда проходит через начальную точку B0[b, f(b)]. Для получе$
ния формулы, определяющей последующие приближения, рассмот$
рим переход от xn к xn+1. В этом случае уравнение хорды BnB0 как
прямой, проходящей через точки Bn, B0, имеет вид
y 1 f (xn )
x 1 xn
2
f (b) 1 f (xn ) b 1 xn .
15
Если для определения xn+1 положить y(xn+1) = 0, то получим
b 3 xn
(n = 0, 1, 2,...).
(1.17)
f 1 b 2 3 f 1 xn 2
Согласно требованиям «а», «б», «в» (см. п. 1.1), наложенным на
функцию f(x) для оценки погрешностей вычислений, используется
неравенство
xn 11 4 xn 3 f 1 xn 2
xn 21 1 x1 2
M 1m
xn 21 1 xn ,
m
где 0< m £ |f'(x)| £ M < 1.
Если при этом M £ 2m, то | xn+1–x*| £ | xn+1–xn| и для заданной
погрешности e вычисления прекращаются при | xn+1 – xn| £ e (как это
имело место и для методов последовательных приближений и метода
касательных).
a)
в)
б)
г)
Рис. 1.9. Метод хорд
16
При выполнении упомянутых требований («а», «б», «в») возмож$
ны и иные картины построений для метода хорд, определяемые соче$
таниями знаков производных f'(x) и f"(x).
Рис. 1.9, а соответствует рассмотренному уже случаю f'(x)>0,
f"(x)>0 (функция f(x) монотонно возрастает и выпукла вниз). Слу$
чай f'(x), f"(x)< 0 (рис. 1.9, б) приводит к аналогичным построени$
ям, и последовательность x1, x2,... оказывается так же возрастаю$
щей.
Однако в случаях f'(x)>0, f"(x)< 0 (рис. 1.9, в) и f'(x) < 0, f"(x)>0
(рис. 1.9, г) после определения каждого xn различными оказываются
знаки значений функций f(a) и f(xn) (а не f(xn) и f(b), как ранее). По$
этому «неподвижной» для всех хорд оказывается точка A0[a, f(a)] (а
не B0[b, f(b)]). В результате расчетными являются формулы
xn 3 a
(n = 0, 1, 2,...),
(1.18)
f 1 xn 2 3 f 1 a 2
а последовательность x1, x2,... оказывается убывающей.
Таким образом, если f'(x)f"(x)>0, то следует использовать форму$
лы (1.17), выбирая за начальное значение x0 = a, если же f'(x)f"(x)< 0,
то используются формулы (1.18) и начальным является x0 = b.
Пример 1.6
Для уравнения f(x) = x3 – 0.2x2 – 0.2x – 1.2 = 0 найти корень x*,
лежащий в интервале [1, 1.5] с точностью до 0.002.
Решение. Так как по условию задачи в указанном интервале ле$
жит один корень x*, проведем его уточнение, пользуясь методом хорд:
f'(x) = 3x2 – 0.4x0.2; f"(x) = 6x – 0.4;
xn 11 4 xn 3 f 1 xn 2
f(1) = – 0.6 < 0; f(1.5) = 1.425 > 0
при
1.0 £ x £ 1.5,
2.4 £ f'(x) £ 5.95; 5.6 £ f"(x) £ 8.6.
Следовательно, f'(x)×f"(x)> 0 при 1.0 £ x £ 1.5.
Пользуясь формулами (1.17) и принимая x0 = a = 1, последова$
тельно находим (b = 1.5):
x1 4
0.6 11.5 3 12
1.425 5 0.6
x2 2 1.15 3 0.173
5 1 4 1.15 ; f(x1) = – 0.173;
1.5 1 1.15
2 1.190 ; f(x2) = – 0.036;
1.425 3 0.173
17
x3 2 1.190 3 0.036
x4 2 1.198 3 0.0172
1.5 1 1.190
2 1.198 ; f(x3) = – 0.0172;
1.425 3 0.036
1.5 1 1.198
2 1.199 ; f(x4) = – 0.0061.
1.425 3 0.0172
Так как |x4 – x3| = |1.199 – 1.198| £ 0.002, то за приближенное
значение корня следует принять x»*x4 = 1.199. (Точное значение кор$
ня x* = 1.2).
1.8. Модифицированный метод хорд
При выполнении требований «а», «б», «в» (см. п. 1.1) относи$
тельно функции f(x) в некоторых случаях быстрее приводит к ре$
зультату (требует меньшего числа последовательных приближений)
модифицированный метод хорд, отличающийся от предыдущего тем,
что каждая новая хорда проводит$
ся не через точки B0 и Bn (или A0 и
An), а через точку Bn[xn, f(xn)] и
точку Bn–1[xn–1, f(xn–1)], отвечаю$
щие предыдущему приближению
(рис. 1.10).
Такой метод, как видно, оказы$
вается близким к методу каса$
тельных, если касательная, про$
водимая в точке Bn, заменяется
хордой, проходящей через эту точ$
ку и предыдущую Bn–1.
Рис. 1.10. Модифицированный
Соответствующие формулы ме$
метод хорд
тода получаются из формул (1.11)
метода касательных при замене значения производной f'(xn) ее при$
ближенным значением
f 4 1 xn 2 5
f 1 xn 2 3 f 1 xn 11 2
xn 3 xn 11
.
Результат получается следующим:
xn 21 4 xn 3 f 1 xn 2
xn 3 xn 11
f 1 xn 2 3 f 1 xn 11 2 (n = 1, 2, …).
(1.19)
Условие окончания вычислений для достижения заданной точно$
сти имеет тот же вид (1.14) | xn+1 – xn| £ e.
18
Правило выбора начального приближения сохраняется тем же,
что и в методе касательных. А именно, в качестве начальной точки
выбирается тот из концов интервала [a, b], в котором f(x)f"(x) > 0.
Для проведения хорды теперь необходимо задать еще одно значение
x1, расположенное внутри интервала [a, b] так, чтобы имело место и
f(x0)f(x1) > 0.
1.9. Комбинированный метод (вариант 1)
Данный вариант метода предполагает объединение метода каса$
тельных и метода хорд.
Пусть на интервале [a, b] f'(x) > 0, f"(x) > 0 (рис. 1.11). Тогда
применение метода касательных (рис. 1.7) дает убывающую после$
довательность приближений x1, x2 ,1 , сходящуюся к искомому кор$
ню x*. В то же время применение метода хорд (рис. 1.9, а) дает возра$
стающую последовательность x1, x2,..., имеющую тот же предел x*.
Следовательно, объединение этих методов, состоящее в их последо$
вательном одновременном применении, дает возможность параллель$
ного вычисления корня x* с избытком и с недостатком. В частности,
значащие цифры, общие для xn и xn, принадлежат и точному значе$
нию корня x*.
Для рассматриваемого случая выберем для метода касательных
начальное значение x0 1 b и для метода хорд –x0 = a. Тогда последо$
вательные приближения по методу касательных определяются фор$
мулами (1.11)
xn 11 3 xn 4
f 1 xn 2
f 5 1 xn 2 (n = 0, 1, 2, …),
(1.20)
а по методу хорд – формулами (1.17):
xn 11 4 xn 3 f 1 xn 2
b 3 xn
f 1 b 2 3 f 1 xn 2
(n = 0, 1, 2, …).
На каждом шаге вычислений погрешность может быть уменьше$
на, если принять
1
1 xn 4 xn 2 .
2
Следовательно, процесс вычислений может быть прекращен,
если xn 11 1 xn 11 2 3 , где e – заданная точность вычислений.
Описанная процедура графически показана на рис. 1.11.
Рассмотренный метод относится к случаю f'(x)>0, f"(x)>0. При
других вариантах этих неравенств для метода касательных различ$
x1 3
19
ным образом выбирается x0, а для
метода хорд – как x0, так и сами фор$
мулы (1.17) и (1.18).
В результате приходим к следу$
ющим утверждениям.
При использовании данного ком$
бинированного метода необходимо
поступать следующим образом.
В методе касательных исполь$
зуются единые формулы (1.11) и
за x0 выбирается то граничное зна$
Рис. 1.11. Комбинированный
чение a или b, для которого
метод (вариант 1)
f 3 1 x0 2 4 f 33 1 x0 2 1 0 .
В методе хорд, если f'(x)f"(x)>0, то используются формулы (1.17)
при начальном значении x0 = b; если же f'(x)f"(x)< 0, то формулы
(1.18) при начальном значении x0 = a.
Пример 1.7
Вычислить наименьший положительный корень уравнения ex –
3x = 0 с точностью e = 0.0001, пользуясь комбинированным методом
(вариант 1).
Р е ш е н и е . Как показано в примере 1.5 (п. 1.5), при решении
того же уравнения методом касательных получена возрастающая пос$
ледовательность приближений (табл. 1.1):
x0 1 0; x1 1 0.5; x2 1 0.61; x3 1 0.6190; x4 1 0.61909.
Применяя теперь метод хорд, пользуясь формулами (1.17) и вы$
бирая начальную точку x0 = 1, вычисляем
x1 = 0.78021;
x4 = 0.62398;
x7 = 0.61918;
x2 = 0.67335;
x5 = 0.62050;
x8 = 0.61909;
x3 = 0.63568;
x6 = 0.61948;
x9 = 0.61906.
Как видно, для x4 и x8 совпадают 5 значащих цифр. Следова$
тельно, в качестве искомого корня с требуемой точностью можно при$
нять x*= 0, 61909.
1.10. Комбинированный метод (вариант 2)
Этот вариант метода предполагает объединение метода касатель$
ных и модифицированного метода хорд. Его особенность, по сравне$
нию с предыдущим вариантом, заключается в том, что на каждом
шаге метод хорд применяется к новому интервалу 1 xn , xn 2 (см. п. 1.8).
20
При этом сохраняется без изме$
нений все, относящееся к методу ка$
сательных и построению последо$
вательных приближений x0, x1,...
(см. п. 1.3). В то же время каждая
последующая хорда, пересечение
которой с осью абсцисс определяет
приближение xn+1, проводится че$
рез точки с абсциссами xn и xn, от$
вечающие предыдущему приближе$
нию (рис. 1.12).
Соответствующие формулы ме$
тода принимают вместо (1.19) сле$
дующий вид:
xn 11 4 xn 3 f 1 xn 2
Рис. 1.12. Комбинированный
метод (вариант 2)
xn 3 xn
f 1 xn 2 3 f 1 xn 2 (n = 1, 2, …).
(1.20)
Выбор начальных точек для каждого из составляющих методов
производится так, как это было указано в пп. 1.5 и 1.8 соответствен$
но.
Пример 1.8
Вычислить с точностью e = 0.0001 положительный корень урав$
нения f(x) = x2 – sin 5x = 0.
Р е ш е н и е . Для определения искомого корня представим урав$
нение в виде x2 = sin 5x и построим график функций y = x2, y = sin 5x
(рис. 1.13).
Как видно, существует единственный положительный корень x*,
который можно считать заключенным в интервале [0.5;0.6].
Для уточнения корня воспользуемся описанным комбинированным
методом, полагая a = 0.5, b = 0.6.
Так как f'(x) = 2x – 5cos5x,
f"(x) = 2+25sin5x и f(b)f"(b) = (0.62 –
– sin(5×0.6))(2+25sin(5×0.6)) > 0, то
для нахождения последовательных
приближений x1, x2, 1 методом ка$
сательных выбираем согласно (п. 1.5),
x0 = b = 0.6. Сами приближения вы$
числяем, пользуясь формулами
(1.20).
Так как имеет место и f(x) f"(x) > 0
Рис. 1.13. Отделение корней
уравнения
при 0.5£ x £ 0.6, то согласно (п.1.7),
21
та же точка x0 = 0.6 является начальной для построения последова$
тельности приближений x1, x2,… модифицированным методом хорд.
Эти приближения находим, пользуясь формулами (1.20). Результа$
ты вычислений сведены в табл. 1.2.
Таблица 1.2
Приложения
x
x2
sin 5x
f(x)
f'(x)
a
0.5
0.25
0.59847
–0.34847
–
b
0.6
0.36
0.14112
0.21888
6.14995
x0
0.56441
0.31856
0.31413
0.00443
5.87572
x1
0.56142
0.31519
0.32829
–0.01310
–
x2
0.56366
0.31771
0.31769
0.00002
–
x2
0.56360
0.31765
0.31797
–0.00032
–
Как видно из табл. 1.2, значение искомого корня можно принять
равным x* » x = 0.5636, так как при у, x2 и x2 совпадает необходимое
число значащих цифр и погрешность x2 1 x2 2 0.00006 1 3 2 0.0001 не
превышает заданной.
1.11. Комбинированный метод (вариант 3)
В этом варианте предполагается объединение упрощенного мето$
да касательных (см. п. 1.6) и метода хорд (см. п. 1.8). Отличие от
предыдущего варианта заключается в использовании на каждом шаге
значений производной f'(x) (в методе касательных), вычисляемой не
на каждой точке xn, а в одной и той
же начальной точке x0 (рис. 1.14).
Расчетные формулы принима$
ют вид (1.16) для метода касатель$
ных и формулу (1.19) – для мето$
да хорд. Все относящееся к выбо$
ру начальных точек сохраняется
таким же, что и для каждого из
указанных методов.
Пример решения задачи этим
методом полностью совпадает с
предыдущим примером 1.8 (см.
Рис. 1.14. Комбинированный метод
п. 1.10) с той лишь разницей, что
(вариант 3)
22
при вычислении производных (последний столбец табл. 1.2) значе$
ния производной f'(x) для всех точек принимаются равными началь$
ному f'(x) » f'(b) = 6.14995.
1.12. Метод последовательных приближений
для системы двух нелинейных уравнений
Пусть имеется система двух уравнений с двумя неизвестными
53 f 1 x, y 2 4 0,
6
(1.21)
57 g 1 x, y 2 4 0.
Ее решением называется такая пара чисел (x*, y*), которая обра$
щает уравнения (1.21) в тождество: f(x*, y*) = g(x*, y*) = 0. Каждое из
уравнений (1.21) можно «частично» решить относительно x или от$
носительно y, т. е. представить систему в виде:
63 x 4 5 1 x, y 2,
7
69 y 4 8 1 x, y 2.
К такой системе может быть применен метод последовательных
приближений, описанный в п. 1.3.
Выбрав начальное приближение (x0, y0), строим последователь$
ность числовых пар:
x1 = j (x0, y0),
x2 = j (x1, y1),
xn = j (xn–1, yn–1);
y1 = y (x0, y0),
y2 = y (x1, y1),
yn = y (xn–1, yn–1).
При выполнении определенных требований к функциям j и y обе
последовательности xn и yn сходятся: xn® x*, yn® y*.
Можно доказать, что в этом случае пара (x*, y*) будет решением
исходной системы.
Условия, налагаемые на функции j и y и обеспечивающие сходи$
мость, здесь не рассматриваются. Все предлагаемые для решения зада$
чи подобраны таким образом, что последовательности приближений
сходятся. Более того, в этих задачах уравнения системы решены отно$
сительно x и y не «частично», а полностью, т. е. система имеет вид:
63 x 4 5 1 y 2,
7
69 y 4 8 1 x 2.
Процесс решения приобретает в этом случае попеременный харак$
тер. x® y® x® y®…: или y® x® y® x®…
Например, если начать с некоторого значения x0, то последова$
тельность примет вид
23
y1 = y(x0), x1 = j(y1), y2 = y(x1), x2 = j(y2), …
(1.22)
Пример 1.9
Найти методом последовательных приближений корни системы
2
3
4 y 5
66x 7 1.63 1 8 7 9 1 y 2,
2.08 6
6 y 7 ln x 7 1 x 2.
(1.23)
Построив графики обеих функций, видим, что в качестве началь$
ных приближений можно взять x0 = 1.
Рис. 1.15. Предварительная оценка решения системы (1.23)
Уточнение этих приближений по формулам (1.22) дает значения
x* = 1.679927, y* = 0.518751.
1.13. Методические указания к составлению алгоритма реше%
ния нелинейного уравнения (системы уравнений)
Процесс отыскания значения корня уравнения (или корней систе$
мы уравнений) с требуемой точностью сводится к последовательному
вычислению значений рекуррентной числовой последовательности,
каждый i$й элемент которой является i$м приближением к истинно$
му значению корня.
Рекуррентной называется такая последовательность, каждый
член которой получается из предыдущего по одной и той же формуле
xi+1 = f(xi), где f – некоторая функция.
В большинстве рассматриваемых ниже задач последовательности
xi при i® ¥ сходятся к истинному значению корня.
Вычисления могут быть прекращены в том случае, когда два со$
седних члена последовательности достаточно близки: |xi+1 – xi| < e,
где e – заранее заданное положительное число.
24
В нескольких задачах, решаемых методом последовательных при$
ближений с использованием обратных функций, заранее не извест$
но, какая из двух числовых последовательностей сходится, а какая
– расходится. В этом случае необходимо каждый раз осуществлять
две проверки: на сходимость (т. е. достаточную близость соседних
членов) и на отсутствие сходимости (например, рост абсолютной ве$
личины разности соседних членов).
1.14. Варианты заданий
для выполнения лабораторной работы
Таблица 1.3
№
Уравнение
Метод численного
решения, точность
Параметры
1
a
1 cx 2 0
b 3 x3
Дихотомии e = 10–3
a = 1.05; b = 0.1;
c = 2.03
2
ax+bsinx+c = 0
Итераций e = 10–4
a = 4.02; b = 0.2;
c = –0.33
3
tg(ax+b)+cx2 = 0
Ньютона e = 10–3
a = 3.01; b = 4; c = –1
4
ax+bxsinx = 0
Хорд e = 10–5
a = 2.01; b = –1
5
asin2x+blgx = 0
Kомбинированный
e = 3·10–4
a = 2.43; b = –0.573
6
a 1 x 2 2 bx2 3 0
Дихотомии e = 4·10–5
a = 1.23; b = –3.14
7
a+bx2+cex = 0
Итераций e = 5·10$–3
a = 1.03; b = 1.01;
c = –2.33
8
sinax+ bx3+cx = 0
Ньютона e = 6·10–4
a = 2.23; b = –3.14;
c = 1.02
9
ax3 1 b 1 c x 1 2 2 0
Хорд e = 7·10–5
a = 1.11; b = –10.11;
c = –2.02
10
aln(b+cx)+dx3 = 0
Kомбинированный
e = 8·10–3
a = 1.23; b = 2.14;
c = 1.001; d = 0.53
11
2x1a3 x
6 b 7 4c 1 d 5 0
8
9
Дихотомии e = 9·10–4
a = 0.1; b = 2.23;
c = 2; d = –1.03
12
ax+b+cx = 0
Итераций e = 10–4
a = 3.51; b = 1.47;
c = 2.004
13
(a+x)2+bx+c = 0
Ньютона e = 10–5
a = –2.13; b = 1.47;
c = –3.14
3
25
Окончание табл. 1.3
№
Уравнение
Метод численного
решения, точность
Параметры
14
acos(x+b)+cx3 = 0
Хорд e = 2·10–4
a = 2.13; b = 3.62;
c = –4.12
15
tgx2+ax+b = 0
Kомбинированный
e = 2·10–5
a = 3.01; b = 2.017
16
(x+a)5+bx = 0
Дихотомии e = 3·10–4
a = 0.29; b = 2
17
ax2 1 b 1 c x 1 d 2 0
Итераций e = 10–3
a = 1.02; b = –4.37;
c = –2.1; d = –0.5
18
xa (bx+c)2 – 14 = 0
Ньютона e = 4·10–4
a = 3.23; b = 1.2;
c = 3.22
19
ln(x+a)+(x+b)5 = 0
Хорд e = 3·10–5
a = 2.11; b = –4.03
20
ax – bx+c = 0
Kомбинированный
e = 2·10–3
a = 2.03; b = 3.02;
c = 2.73
21
(x+a)2 –ebx = 0
Дихотомии e = 10–4
a = –0.4; b = 0.53
22
(a–bsinx)2+clnx = 0
Итераций
e = 0,5·10–3
a = 1.23; b = 1.45;
c = 1.004
23
(x+a)3+bsincx = 0
Ньютона e = 2·10–4
a = –3.21; b = –1.45;
c = 2.12
24
ax2cosbx – cx = 0
Хорд e = 7·10–5
a = 2.93; b = 3.01;
c = 2.1
25
ax2sinbx– lg(x+c) = 0
Kомбинированный
e = 0,5·10–4
a = 1.003; b = 1.008;
c = 4.15
Дихотомии
e = 2·10–5
a = 2.07; b = 1.19;
c = 1.13; d = –1.001
26
a
bx 1 c
1 d cos bx 2 0
27
ax
1 d lg(ax) 2 0
b 1 cx2
Итераций
e = 10–5
a = 1.01; b = 0.98;
c = 2.03; d = –2.04
28
a sin x 1 b lg x 2 0
Ньютона
e = 4·10–5
a = 2.06; b = –1.06
Хорд e = 5·10–4
a = 2.37; b = –0.99;
c = 0.56
Kомбинированный
e = 6·10–4
a = 4.01; b = 0.97;
c = 2.98; d = –2.98
29
30
26
a
1 becx 2 0
x
acos(bx+c)+dx = 0
Лабораторная работа№ 2
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ
2.1. Постановка задачи
Пусть требуется найти определенный интеграл
b
3 f 1 x 2 dx,
(2.1)
a
где функция f(x) непрерывна на отрезке [a, b]. Если функция f(x)
задана в явном виде и можно найти неопределенный интеграл F(x),
то определенный интеграл вычисляется по формуле
b
3a f (x)dx 1 F(b) 2 F(a)
(2.2)
(формула Ньютона$Лейбница).
Если же неопределенный интеграл данной функции найти слож$
но или функция f(x) задана графически или таблицей, то для вычис$
ления определенного интеграла применяют приближенные форму$
лы. Для приближенного вычисления интеграла (2.1) существует
много численных методов, из которых выделим три:
1) метод прямоугольников;
2) метод трапеций;
3) метод Симпсона.
При вычислении интеграла следует 1
помнить, каков геометрический смысл
определенного интеграла. Если f(x)³ 0 на
отрезке [a, b], то
b
1a f (x)dx
2
3
4
численно равен площади фигуры, ограни$
Рис.2.1. Геометрический
ченной графиком функции y = f(x), от$ смысл численного интеграла
резком оси абсцисс, прямой x = a и пря$
мой x = b (рис. 2.1).
Таким образом, вычисление интеграла равносильно вычислению
площади криволинейной трапеции.
2.2. Метод прямоугольников
Разделим отрезок [a, b] на n равных частей, т. е. на n элементар$
ных отрезков. Длина каждого элементарного отрезка h = (b–a)/n.
27
Точки деления будут: x0 = a; x1 = a+h; x2 = a+2h, …, xn–1 = a+(n–1)h;
xn = b. Эти числа будем называть узлами. Вычислим значения функ$
ции f(x) в узлах, обозначим их y0, y1, y2, …, yn. Стало быть, y0 = f(a),
y1 = f(x1), …, yn = f(b). Числа y0, y1, y2, …, yn суть ординаты точек графи$
ка функции, соответствующих абсциссам x0, x1, x2, …, xn (рис. 2.2). Из
рис. 2.2 следует, что площадь криволинейной трапеции приближен$
но заменяется площадью многоугольника, составленного из n пря$
моугольников. Таким образом, вычисление определенного интеграла
сводится к нахождению суммы n элементарных прямоугольников.
51
4121
52
42
12222
41
51
41
3 3212
11
1121 2 3211
1
Рис. 2.2. Геометрическая иллюстрация метода прямоугольников
b
S 1 3 f (x)dx 1 y0 * h 2 y1 * h 2 y2 * h 2 y3 * h 2 1 2 yn * h 1
a
1 h * (y0 2 y1 2 y2 2 y3 2 1 2 yn 11 );
(2.3)
b
S 1 3 f (x)dx 1 y1 * h 2 y2 * h 2 y3 * h 2 y4 * h 2 1 2 yn * h 1
a
1 h * ( y1 2 y2 2 y3 2 y4 2 1 2 yn ).
(2.4)
Формула (2.3) называется формулой левых прямоугольников, а
(2.4) – формулой правых прямоугольников.
h2
1
S 3 h 9 y 5 x1 4 6 –
28
7
формула средних прямоугольников.
(2.5)
2.3. Формула трапеций
Формула трапеций имеет вид:
b
S 1 3 f (x)dx 1 h * ((y0 2 yn )/2 2 y1 2 y2 2 1 2 yn 11 ).
a
28
(2.6)
Формула (2.6) означает, что площадь криволинейной трапеции
заменяется площадью многоугольника, составленного из n трапеций
(рис. 2.3); при этом кривая заменяется вписанной в нее ломаной.
1
4
51
1
53
52
43 42
13
3 1213 12
41
51
2 1211
11
1
Рис. 2.3. Геометрическая иллюстрация метода трапеции
2.4. Метод Симпсона
Геометрически иллюстрация формулы Симпсона состоит в том,
что на каждом из сдвоенных частичных отрезков заменяем дугу дан$
ной кривой дугой графика квадратного трехчлена.
Разобьем отрезок интегрирования [a, b] на 2n равных частей дли$
ны h = (b–a)/2n. Обозначим точки разбиения x0 = a; x1 = x0+h, …, xi =
x0+i*h, …, x2n = b. Значения функции f в точках xi обозначим yi, т. е.
yi = f(xi). Тогда согласно методу Симпсона:
b
S 3 f 1 x 2 dx 3 1 b 4 a 2 /6n 5 1 y0 6 4y1 6 2y2 6 1 6 4y2n 11 6 y2n 2 7
a
2n 11
8
9
i 11
7 1 b 4 a 2 /6n 5 y0 6 y2n 6 3 6 1 412
5 yi .
i 20
2
1
(2.7)
2.5. Задания для выполнения лабораторной работы
0.8
1.
2
0.2
1.3
2.
2
sin(2x 1 0.5)dx
2 1 cos(x2 1 1)
sin(0.5x 1 0.4)dx
0.5 1.2 1 cos(x
2
1 0.4)
1.4
7.
x2 1 5dx
2 2.3 1 sin(1.5x 1 1)
0.6
2.1
8.
3
1.3
sin(x2 1 1)dx
22 x
29
0.9
3.
2
0.3 1.5 1 sin(x
1.0
4.
2
2
2
0.6
1 0.6)
sin(x 1 1.4)dx
0.4 0.8 1 cos(2x
1.0
5.
cos(0.8x 1 1.2)dx
2
1 0.5)
cos(0.6x2 1 0.4)dx
1.4 1 sin2 (x 1 0.7)
1.2
tg (x2 )dx
9. 2
x 11
0.5
0.63
30
cos(0.4x2 1 0.5)dx
2 sin(x 1 0.7) 1 x
0.6
x 1 1lg(x 1 3)dx
0.15
0.8
11.
2
0.4
1.0
6.
2
10.
1.6
12.
2
0.8
tg (x2 1 0.5)dx
1 1 2x2
0.3x2 1 2.3dx
1.8 1 2x 1 1.6
Лабораторная работа № 3
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
3.1. Общие сведения о дифференциальных уравнениях
Обыкновенное дифференциальное уравнение (ОДУ) имеет общий
вид
Ф(x, y, y', y",..., y(n)) = 0,
где x – независимая переменная; y – искомая неизвестная функция
этой переменной; y', y" – производные у по x. Наивысший порядок n
этих производных определяет порядок ОДУ. Здесь будут рассматри$
ваться ОДУ первого порядка в виде
y' = f(x, y).
(3.1)
Решением ОДУ (3.1) называется функция y = y(x), обращающая
его в тождество. Геометрически на плоскости ОXY решению ОДУ со$
ответствует семейство интегральных кривых y = j(x, C) (где C – про$
извольная постоянная), определяющее общее решение ОДУ. Част$
ные решения получаются из общего выбором конкретного значения
постоянной C, при этом из семейства интегральных кривых выделя$
ется одна кривая.
Нахождение частного решения уравнения (3.1), удовлетворяю$
щего начальному условию
y(x0) = y0,
(3.2)
называют задачей Коши, которая сводится к определению интеграль$
ной кривой y = y(x), проходящей через заданную начальную точку
М0 (рис. 3.1, а).
Найти аналитическое решение ОДУ удается не всегда; для мно$
гих случаев такое решение невозможно. Поэтому практически важ$
ны численные методы решения ОДУ. При численном решении за$
данный интервал независимой переменной [x0, xk] разбивается на
равные промежутки шириной h, т. е. выбирается система точек xi
= x0+ih (где i = 0, 1, 2,..., k). Этим дискретным значениям x i соот$
ветствуют дискретные значения yi искомой функции y(x) (рис. 3.1,
а). В каждой точке Mi(xi, yi) наклон касательной должен соответ$
ствовать заданному ОДУ (3.1), т. е. y'i = f(x i, yi). В результате по$
лучается решение в виде таблицы (рис. 3.1, б), т. е. в отличие от
31
Рис. 3.1 Численное решение ОДУ: a – график решения; б – таблица
аналитического решения (кривой y = y(x)) численные методы дают
всегда решение в дискретной форме.
На каждом шаге решения требуется переход от точки i к точке
i+1. Для независимой переменной очевидно, что xi+1 = xi+h. Точность
полученного решения зависит как от выбранного метода интегриро$
вания, так и от величины шага h. C уменьшением шага h точность
решения возрастает, однако возрастает и объем вычислений. На прак$
тике нередко сравнивают решение, полученное при шаге h, с решени$
ем, полученным при шаге h/2. Если значения полученных решений
совпадают с заданной точностью, то задача численного решения ОДУ
считается выполненной. Вопросы оценки точности решения при чис$
ленных методах имеют первостепенное значение.
3.2. Классификация численных методов решения
обыкновенных дифференциальных уравнений
Будем исходить из уравнения
Y = F(Y, X, L, t).
(3.3)
Для простоты рассмотрим сначала скалярное уравнение
ў = f(y, t).
(3.4)
Численным решением уравнения является табл. 3.1, которая на
каждом шаге h дает численные значения переменной ў(kh), достаточ$
но близкие к истинным значениям y(kh).
Классические методы численного интегрирования дифференциаль$
ных уравнений восходят к Эйлеру (XVII в.), Адамсу (вторая полови$
на XIX в.). Появление ЭВМ и бурное их развитие для решения широ$
32
Таблица 3.1
k
t
y
0
t0
ў0
1
t0 + h
ў1
2
t0 + 2h
ў2
…
…
…
k…
t0 + kh…
ўk
кого круга задач науки и техники потребовало переоценки класси$
ческих методов численного интегрирования, их машинной ориента$
ции и разработки новых, собственно машинных методов, учитываю$
щих специфику ЭВМ. Машинная ориентация классических числен$
ных методов состоит в распространении этих методов на системы урав$
нений большой размерности, автоматизации операций реализации
методов на ЭВМ, выборе наиболее подходящего метода для решения
конкретной задачи и величины шага в этом методе.
Аналитикочисленные методы, основанные на классическом раз$
ложении правых частей в степенные ряды Тэйлора в аналитической
форме, обычно исключаются из практического применения ввиду
необходимости вычислений производных функций порядка n>1. Од$
нако применение средств аналитических преобразований на ЭВМ по$
зволяет для ряда ММ «заготавливать» эти производные в аналити$
ческом виде, что принципиально изменяет численные методы на сте$
пенных разложениях, делает их эффективными. В соответствии с
отмеченными обстоятельствами рассмотрим построение таких мето$
дов на основе разложения Тэйлора.
Будем искать решения y(t) в форме ряда
1
Y 1 3 ak tk 2 a0 , a0 1 Y0 ,
(3.5)
k 21
сходящегося на промежутке 0 £ t £ R. Тогда все дело сводится к
нахождению неопределенных коэффициентов ak, k Î [0, ¥]. Для их
нахождения подставим (3.5) в (3.4), предварительно разложив f(y,
t) в ряд Тэйлора относительно точки t = 0, y = y0,
1
3f
3f
1 32 f
7 kaktk21 4 f 1 0,y0 2 5 3t t 5 3t 1 y 6 y0 2 5 2! 3t2 t2 5
k 31
5
1 32 f
1 32 f
2
y
6
y
5
1
2
1 y 6 y0 2 t 5 1
0
2! 3y2
2! 3y3t
(3.5')
33
Приравнивая в этом разложении коэффициенты при одинаковых
степенях t, найдем
a0 3 y0 ,
a1 3 f 1 0, y0 2,
a2 3
1 4 6f 6f 5
8 7 a1 9
2! 6t 6y y
,
0 ,0
a3 3
5
1 4 1 62 f 6f
1 62 f 2 1 62 f
7 a2 7
a1 7
a1 9
88
2
2
6y
3! 2! 6t
2! 6y
2! 6y6t 9
,
y0 ,0
1111111111111111111 .
(3.6)
Таким образом, можно последовательно найти ak и тем самым по$
строить решение y(t).
Однако этот алгоритм обладает существенными недостатками –
необходимостью ограничиваться только промежутком сходимости
ряда (3.5) – R, оперировать с большим числом членов этого ряда и
брать производные высокого порядка от правой части в (3.5'). От
первого недостатка можно избавиться, воспользовавшись аналити$
ческим продолжением решения y(t)
1
k
y 3 6 ak 1 t 4 tj 2 5 y0 ,
k 21
где tj – переменный центр разложения, а(t – tj) < R0, j = 1, 2, 3, …, tn = 0.
В этом случае коэффициенты разложения ak каждый раз пересчи$
тываются для нового центра разложения tj, yi (рис. 3.2, а).
От второго недостатка избавляются применением асимптотичес$
кого представления y(t),
n
3 y 1 t 2 3 6 ak 1 t 3 tj 2 4 5,
k
(3.6')
k 11
Для асимптотического представления y(t) удобно воспользовать$
ся формулой Тэйлора в точке t = tj:
n
y 1 t 2 6 y 1 tj 2 7 y
k 21
1 i2
1 tj 2
1 t 3 tj 2
i!
i
7 O 41 t 3 tj 2
8
n 11 5
,
9
где O[(t – tj) n+1] – остаточный член, представляющий малую по (t – tj)
выше n$го порядка y(i)(tj), i$я производная от y(tj)|t = ti. В теории рядов
34
Тэйлора доказано, что если для какой$либо функции j(t), имеющей
производные до n$го порядка включительно, в точке t = h выполня$
ется условие
11 1 h 2 4 2 4 31 n 2 1 h 2 4 0 ,
3 1 h 2 4 31 1 h 2 4 3
то j(t) ® 0 при t® h и представляет собой бесконечно малую выше
n$го порядка малости. Тогда, обозначив
1i 2
n
y
n
hi
31 2 1 h 2 4 y 5 y j 5 6
i
!
i 11
и обеспечив условие (3.6') выбором n и h, получим
1
1
1
2
3 4 O hn 11 5 yn 11 1 6h 2
R n 11
,
1 n 7 12 !
2
1
1
323
31
1
1124
3
42115313678794263
11
1
31
1124 761194263
11
1124
1124
4
3164561194
64367879126
11
1
3124
64 1
3124
31
71124 761191263
Рис. 3.2. Операции построения численных методов РунгеIКутты
35
где Q < 1.
На промежутке TО (tj+1, tj) асимптотическое разложение y(t) при$
мет вид
n
y 1 tj 11 2 3 yi 1 tj 2 4 5 y
1i2
i 21
i
1 tj 2 hi !
.
В этом выражении есть два регулируемых параметра n и h. Выбо$
ром этих параметров можно добиваться точности порядка h при чис$
ле членов разложения n. Современные средства программирования
позволяют автоматически получать аналитические выражения для
ai, в (3.6). Поэтому процедуру численного интегрирования можно по$
строить в соответствии с алгоритмом (3.6), возложив на ЭВМ опреде$
ление производных, расчет коэффициентов ai, выбор шага h ´R. Рас$
пространение аналитико$численного алгоритма (3.6), составленно$
го для скалярного управления на систему уравнений, состоит в вы$
числении производных и составлении в аналитическом виде коэффи$
циентов для векторного уравнения
1
2
Y j 11 6 Y j 7 hF Y j ,tj 7
h 2 4 3F 5
3F
* Y1 7
8
9
3t
2 3Y Y ,t
j
j
* Y1 7 2,
Yj ,tj
где ¶F/¶t – вектор, полученный покомпонентным частным диффе$
ренцированием, а ¶F/¶ Y есть (n1 ´ n1) матрица Якоби функции F(Y,
t) системы ( 3.3 ).
Построение численных методов, в которых переход от ў(t) к ў(tj+1)
не требовал бы вычисления производных от правых частей в (3.4 ),
основано на идее Рунге. Принцип построения таких формул числен$
ного интегрирования можно усмотреть из геометрических построе$
ний (рис. 3.2, б, в). Продолжая этот процесс и организуя комбина$
ции ki(h) для нескольких точек внутри промежутка tj+1 – tj, запишем
общую формулу этой группы численных методов
n
y j 11 3 y j 4 5 bi ki 1 h 2,
i 21
где
ki 1 h 2 5 hf 39tj 6 7 i h, y j 6 8i1k1 1 h 2 6 8i2k2 1 h 2 6 1 6 8i,i 11ki 11 1 h 2 4
. (3.7)
Из формулы (3.7) при различных значениях aj и bj получим как
частные случаи известные формулы численного интегрирования Рун$
ге$Кутты разного порядка точности:
первого порядка –
1
2
y j 11 3 y j 4 k1 1 h 2, k1 1 h 2 3 hf y j ,tj – метод Эйлера;
36
второго порядка –
1
1
y j 11 5 y j 6 37k1 1 h 2 6 k2 1 h 2 48 ; k2 1 h 2 5 hf y j 6 k1 1 h 2 h,tj 6 h
2
Эйлера$Коши;
2
– метод
h
h4
3h4
3h4
3
y j 11 5 y j 6 k2 7 8, k2 7 8 5 hf 7 yj 6 k1 1 h 2,tj 6 8 – усовершенство$
2
2
92
92
9
ванный метод Эйлера;
четвертого порядка –
4
13
3h4
3h4
y j 11 5 y j 6 7 k1 1 h 2 6 2k2 7 8 6 2k3 7 8 6 k4 1 h 2 8,
69
2
2
9 9 3
h 3h4
h4
3h4
k3 7 8 5 hf 7 y j 6 k2 7 8,tj 6 8,
2
2
2
2
9 9 9
3
4
3h4
k4 1 h 2 5 hf 7 y j 6 hk3 7 8,tj 6 h 8.
2
9 9
Точность представления y на одном шаге в каждой формуле оце$
нивается остаточным членом разложения функции j(h) в ряд Тей$
лора:
O 1 31t 22 5
3r 11 1 4h 2
1 r 6 12
!
hr 11 .
Величина r в этом уравнении называется порядком точности ме$
тода, или степенью метода. Важно сразу же отметить, что выбор шага
h, исходя из заданной точности на одном шаге, не гарантирует требу$
емой точности конечного результата.
Следующая группа методов численного интегрирования, осно$
ванная на идеях Адамса, строится на основе аппроксимации подын$
тегральной функции в формуле Ньютона$Лейбница:
yj 11 3 yj 4 5
tj 11
tj
y1 1 t 2 dt –
(3.8)
известного полинома. Эта формула получена после интегрирования
уравнения (3.4) на промежутке tj – tj+1. Если в качестве такого поли$
нома выбирать полином Ньютона, то y(t) можно представить в виде
y1 1 t 2 3 a0 4 a1 1 t 5 tj 2 4 a2 1 t 5 tj 21 t 5 tj 11 2 4 2
(3.9)
37
Задавая предыдущие значения ў'(tj), ў'(tj–1), ў'(tj–2), … (предпола$
гается, что они вычислены заранее), вычисляем неизвестные коэф$
фициенты
a0 4 y1 1 tj 2, a1 4
y1 1 tj 11 2 3 y1 1 tj 2
h
,2,
(3.10)
где h = |tj tj+1|.
Подставляя (3.10) в (3.9), после интегрирования получим
3 y1 1 tj 12 2 5 y1 1 tj 2 4 h 3 y1 1 tj 12 2 5 2y1 1 tj 11 2 4 5h
7 96
7 92 ,
y 1 tj 21 2 8 y 1 tj 2 9 hy1 1 tj 2 9 6
2
12
или в общем виде
1
n
2
y 1 tj 11 2 3 y 1 tj 2 4 h 5 bi y1 1 tj 21 2 4 O hi 11 2 .
i 30
(3.11)
Методы типа (3.11) называют также многошаговыми в отличие
от одношаговых (3.7), поскольку для получения ўj+1 необходимо
предварительно вычислить n предыдущих значений ўj.
Очевидным недостатком методов этого типа является необходи$
мость предварительного вычисления «разгонных» значений ў'(tj+1).
Методы (3.7), (3.11) реализованы в виде стандартного программ$
ного обеспечения современных ЭВМ.
Применимость методов (3.7), (3.11) тесно связана с выбором шага
h в формулах численного интегрирования. В случае одного уравне$
ния (3.4) значение h выбирается из соображений точности на одном
шаге, и оно же обеспечивает устойчивость вычислительного процес$
са. Когда же рассматривается система уравнений (3.3), то именно
устойчивость разностных уравнений (3.7), (3.11) ограничивает шаг.
Наибольшие шаги, обеспечивающие устойчивость разностной схе$
мы, или, как их обычно называют, критические шаги, различают
для разных методов. Так, например, если собственные значения вы$
численной в данной точке t = t* матрицы Якоби системы веществен$
ны, то критический шаг выбирается из соотношения
h1
r
2 max
,
где lmax – максимальное собственное значение; r – число, зависящее
от метода. Для метода Рунге$Кутты четвертой степени оно равно 2.78,
для методов Эйлера и Эйлера$Коши – 2 и для метода Адамса четвер$
той степени – 0.3. Причем даже незначительное увеличение h против
38
допустимого приводит к потере вычислительной устойчивости,
«взрыву погрешностей».
В приведенном далее примере эта особенность численного интег$
рирования прослеживается в явном виде.
Пример 3.1
Численное интегрирование дифференциального уравнения
y1 3 4104 y 5 e 1t , y 1 0 2 3 y0 , 0 6 t 7 tn1n1 3 1 c .
Применим метод Эйлера и положим h1 1
1
tn1n1
n
(3.12)
, тогда
2
yk21 3 yk 4 h 5104 yk 4 e 1 hk ,
1
4 h 1 510 y
2
4 e 1 2,
k 3 0 y1 3 y0 4 h 5104 y0 4 e0 ,
k 3 1 y2 3 y1
4
1
h
1111111111111
1
2
1
2
k 3 n yn 21 3 yn 4 h 5104 yn 4 e 1nh 3 yn 1 5 h104 4 he 1nh .
Подставляя последовательно предыдущее уравнение в последую$
щее, получим
1
yn 11 3 y0 1 4 104 h
2
n 11
1
3
n
2
5 h 6 1 4 104 h e1
i 0
i
n 2i 2 h
.
(3.13)
Решением уравнения (3.12) очевидно является функция
y 1 y0 e 110 t 2 e 1t .
Условие вычислительной устойчивости метода применительно к
(3.13 ) очевидно будет:
4
|1 – 104h| < 1, откуда h<2•10–4.
Следовательно, надо сделать n = 1/2•104 – 5000 шагов для вос$
произведения элементарной функции Y(t) = e–t на отрезке [0, 1]. Не
спасает и нулевое начальное условие Y0 = 0 в (3.13), ибо вынужден$
ная соответствующая в (3.12) также начинает неограниченно расти
при h>2•10–4.
С другой стороны, одной из основных особенностей реальных сис$
тем и устройств является весьма большой количественный разброс
частот колебаний в процессах динамики, что отражается в плохой
39
обусловленности матрицы Якоби соответствующих дифференциаль$
ных уравнений движения. Ввиду этого масштаб времени
машинное время
( реальное время ), определяемый минимальным собственным зна$
чением tреш > | L–1|, оказывается недопустимо большим.
В табл. 3.2 приведены шаг и масштаб времени, полученные в ре$
зультате решения на ЭВМ рассматриваемыми методами тестовой
задачи вида [2]:
5000 X1''+ 30 000 000(X1– X3) = 7000,
14 X2'' – 235 000X3' = 0,
10 X3''+ 235 000X3'+ 30 000 000( X3 – X1) = 0.
X1(0) = X2(0) = X3(0) = X1'(0) = X2'(0) = X3'(0) = 0. Требуемая
точность решения была назначена D = 0, 02, а время одного варианта
tреш = 10 мин. Результаты тестирования численных методов пред$
ставлены в табл. 3.2.
Таблица 3.2
Метод
Шаг интегрирования, с
Масштаб времени
10–5
2•103
Эйлера–Kоши
2•10–5
2•102
Рунге–Kутты
10–4
102
Адамса 4$го порядка
10–5
5•102
Эйлера
3.3. Метод степенного ряда
Если правая часть ОДУ y' = f(x, y) является аналитической функ$
цией, то искомое решение y(x) может быть представлено степенным
рядом
yi 11 4 yi 5 6yi 4 yi 5 hyi3 5
h2
h2 1 n 2
yi33 5 1 5
y i.
2 !
n !
(3.14)
Точность решения (погрешность вычислений) можно оценить ве$
личиной hn отброшенного члена ряда. Но этот метод удобен только в
тех случаях, когда правая часть заданного ОДУ позволяет получить
сравнительно простые выражения для производных высших поряд$
ков. В общем случае этот метод трудно использовать, так как ЭВМ
«не умеет» дифференцировать.
40
Пример 3.2
Найти решение ОДУ y' = x–y при начальном условии y(0) = 1 на
отрезке 0 £ x £ 1 с шагом h = 0.1, используя пять членов ряда (3.14)
(точное решение y = 2e–x+x–1).
Р е ш е н и е . Дифференцируя правую часть заданного ОДУ, получим
y' = x – y; y" = 1– y'; y''' = – y"; yэv = – y'''; yv = – yэv.
Отсюда выражение (3.14) для степенного ряда будет
h2
1 x 5 y 5 12 4
2
h3
h4
h5
4 1 x 5 y 5 12 5
1 x 5 y 5 12 4
1 x 5 y 5 12.
6
24
120
yi 11 3 yi 4 h 1 x 5 y 2 5
Можно убедиться, что такое решение совпадает с точным в пяти
знаках.
3.4. Метод Эйлера и его модификации
Если в степенном ряду (3.14) сохранить только член первого по$
рядка, то получим
yi+1 = yi+hyi' = yi+h·f(xi, yi),
(3.15)
при этом интегральная кривая ОДУ y(x) заменяется ломаной, на$
клон звеньев которой соответствует yi' на каждом участке от yi до
yi+1 (рис. 3.3). Такой «простой» метод Эйлера весьма не точен. Из
рис. 3.3 видно, что погрешность d накапливается, на каждом звене
2
она порядка h . Уменьшение шага h уменьшает погрешность, но не
устраняет ее. Поэтому простой метод Эйлера применяется редко и
обычно для того, чтобы получить первое «грубое» представление о
характере интегральной кривой. Усовершенствованный метод Эйле$
ра$Коши предполагает использование формулы
yi+1(0) = yi+h·f(xi, yi)
(3.16)
(0)
простого метода Эйлера для «предсказания» значения yi+1 (нуле$
вая итерация). Затем производится «уточнение» yi+1(1), при котором
определяется среднее арифметическое значение y' в начале и конце
интервала шириной в h. Следовательно, «уточняющая» формула име$
ет вид (первая итерация)
112
yi 11 5 yi 6
1
2
h3
1k2
f 1 xi , yi 2 6 f xi 11, yi 11 4 .
8
2 79
(3.17)
41
В итерационном методе Эйлера «уточнение» по предыдущей фор$
муле повторяется несколько раз, а именно:
1 k 112
2
1
h3
1k2
f 1 xi , yi 2 6 f xi 11, yi 11 4 ,
(3.18)
8
2 79
где k – номер итерации. Для расчетов на ЭВМ эту формулу удобнее
представить в виде
yi 11
5 yi 6
1 k 112
yi 11
5 Bi 6
1
2
h3
1k2
f xi 11 , yi 11 4 ,
8
2 79
(3.18.a)
h
3 f 1 xi , yi 2 48 не зависит от итераций.
27
Из рис. 3.4 видно, что точка B получается применением простого
метода Эйлера при шаге h/2. Из точки B проводим прямые, соответ$
ствующие касательным в точках нулевой yi+1(0), первой yi+1(1), вто$
рой yi+1(2) (и т. д.) итерациях. Разность последней и предпоследней
итераций d = | yi+1(k+1)– yi+1(k)| не должна превышать допустимой
ошибки e, т. е.
d £ e.
(3.19)
где Bi 5 yi 6
Рис. 3.3. Простой метод Эйлера
Рис. 3.4. Метод Эйлера с итерациями
Если за 3–4 итерации условие (3.19) не достигается, то шаг h сле$
дует уменьшить. Напротив, если заданная точность достигается сра$
зу за одну итерацию (k = 1), то шаг можно увеличить.
Значение e, используемое в качестве критерия окончания итера$
ций, следует отличать от ошибки, с которой получено решение, так
как процесс итераций сходится к значению, которое может не соот$
ветствовать точному решению ОДУ. Метод Эйлера с итерациями до$
пускает ошибку вычислений (точность решения) не хуже h3 (т. е. на
порядок выше, чем простой метод Эйлера).
42
3.5. Метод Рунге – Кутты
Этот метод является одним из наиболее распространенных, так
как точность его сравнительно велика (ошибка порядка h5). Переход
от yi к yi+1 совершается по формуле
h
1T1 4 2T2 4 2T3 4 T4 2 .
6
Здесь используются вспомогательные величины
yi 11 3 yi 4
(3.20)
T1 = f(xi, yi) = y'(M1) = tg a1,
(3.21а)
h
h 4
3
T2 6 f 9 xi 7 , yi 7 T1 6 y5 1 M2 2 6 tg82 ,
2
2 (3.21б)
h
h 4
3
T3 6 f 9 xi 7 , yi 7 T2 6 y5 1 M3 2 6 tg83 ,
2
2 (3.21в)
T4 4 f 1 xi 5 h , yi 5 hT3 2 4 y3 1 M4 2 4 tg64 .
(3.21г)
Геометрический смысл их состоит в том, что они выражают тан$
генс угла наклона касательных, определяемый правой частью ОДУ в
точках М1, М2, М3, М4 внутри интервала шагом h (рис. 3.5, а).
Рис. 3.5. Метод Рунге–Кутты с постоянным (а) и переменным (б) шагами
Автоматический выбор шага. При расчетах на ЭВМ методом Рун$
ге$Кутты нередко пользуются так, что на каждом этапе расчет ведет$
ся двояко: с шагом h и (дважды) с шагом h/2 (рис. 3.5, б). При шаге
h получается результат U, а при шаге h/2 получается результат W.
Модуль разности этих результатов обозначим D = |U–W|. Если D £ e,
то точность достигнута, шаг h не изменяется и следующий этап вы$
полняется аналогично. Если D = e, т. е. точность не достигнута, то
первоначально выбранный шаг уменьшается вдвое и расчет на дан$
43
ном этапе повторяется при этом уменьшенном шаге. Если же ока$
жется, что заданная точность намного (например, в 32 раза) превы$
шает разницу D (т. е. e > 32D), то шаг на следующем этапе можно
увеличить вдвое.
3.6. Метод Адамса
В предыдущих «одношаговых» методах (степенного ряда, Эйле$
ра, Рунге$Кутты) для продвижения решения на один шаг достаточно
было знать данные только в одной пре$
дыдущей точке.
Метод Адамса – многошаговый.
Здесь надо знать результаты в несколь$
ких (обычно четырех) точках xi, xi–1,
xi–2, xi–3, чтобы найти значение реше$
ния ОДУ в следующей точке xi+1, про$
двигаясь вперед на шаг h (рис. 3.6).
В четырех предыдущих точках не$
обходимо знать производные, а имен$
но:
Рис. 3.6. Метод Адамса
A = yi, B = yi–1, C = yi–2, D = yi–3.
Вначале используют «предсказывающую» формулу (нулевая ите$
рация)
102
h
355 A 7 59B 6 37C 7 9D 4
(3.22)
24
(при отсутствии итераций этим и ограничиваются). При использова$
нии итераций далее применяется «уточняющая» формула
yi 11 5 yi 6
1 k 112
1
2
h 3
1k2
19 A 7 5B 6 C 7 9f xi 11, yi 11 4
(3.23)
9
24 8
(k – номер итерации). Если модуль разности результатов вычислений
по формулам (3.22) и (3.23) меньше допустимой ошибки e, то про$
движение на шаг h выполнено и далее переходят к следующему шагу.
Если требуемая точность не достигается, то повторно применяя фор$
мулу (3.23), продолжают процесс итераций. Если за максимально
допустимое число итераций (например, четыре) точность все же не
достигается, то надо уменьшить шаг h вдвое и с этим половинным
шагом (h/2) повторить все решение сначала. Недостаток метода Адам$
са в том, здесь трудно менять шаг «на ходу», так как его надо менять
и в четырех предыдущих точках. Погрешность вычислений по мето$
5
ду Адамса (с итерациями) – величина порядка h (как и в методе Рун$
yi 11
44
5 yi 6
ге$Кутты). Если шаг менять не приходится, то преимуществом мето$
да Адамса является возможно меньшее число обращений к правой
части ОДУ, чем в методе Рунге–Кутта.
3.7. Алгоритмизация методов численного решения ОДУ
Одношаговые методы решения ОДУ
Схемы алгоритмов одношаговых методов решения ОДУ просты и
мало отличаются друг от друга (рис. 3.7). Различие будет состоять
только в блоке 5, где совершается переход от i$й к (i+1)й точке реше$
ния, т. е. yi+1 = yi+Dyi(xi, yi, h), где приращение Dyi является функци$
ей xi, yi, h.
Рис. 3.7. Обобщенная схема алгоритма одношаговых методов решения ОДУ
Для простого метода Эйлера содержимое блока 5 – это не что иное,
как вычисление по формуле (3.15), т. е. Y = Y+H*F(X, Y). В методе
степенного ряда – это вычисление по формуле (3.14), т. е. Y =
Y+H*F(X, Y)+H*H*Y2/2+H**3*Y3/6, (где Y2 и Y3 – арифметичес$
кие выражения для y" и y"' соответственно). Наконец в методе Рун$
ге–Кутты – это вычисление по формуле (3.20) Y =
Y+(T1+2*(T2+T3)+T4)*H/6, где T1, T2, T3, T4 определяются фор$
мулами (3.21).
Метод Рунге – Кутты с переменным шагом
Фрагмент алгоритма этого метода (соответствующий блоку 5 схе$
мы на рис. 3.7) показан на рис. 3.8. Здесь три раза используется об$
ращение к подпрограмме RK(X, Y, H, XN, YN), входными парамет$
рами которой являются входные значения X, Y (в начале интервала
45
Рис. 3.8. Фрагмент схемы алгоритма метода Рунге – Кутты
согласно рис. 3.5, б), а также шаг H. Выходными являются парамет$
ры XN, YN в конце интервала.
Первое обращение к подпрограмме RK при входных параметрах
X, Y, N дает выходные параметры X2, U (см. рис. 3.5, б). Второе
обращение к RK, отличающееся от первого значением шага H/2, по$
рождает промежуточные выходные параметры X1, V в середине ин$
тервала. Третье обращение с выходными параметрами X1, V, H/2
дает выходные параметры X2, W (конец интервала – см. рис. 3.5, б).
Затем вычисляем модуль разности значений U (шаг H) и W (шаг
H/2), т. е. D = |U–W| (блок 4). Если D>EPS (блок 5), то точность еще
не достигнута, и, присвоив U = V, уменьшаем шаг вдвое H = H/2
(блок 8). Снова повторяем операции блоков 1–4 при H = H/2. Если
D £ EPS, то вычисляем окончательно Y = W–(U–W)/15 и присваива$
46
ем X = X2 (блок 6). Если еще при этом EPS ³ 32*D (блок 7), то шаг
можно удвоить (блок 9). Иначе шаг остается неизменным.
Метод Адамса с итерациями
Схема алгоритма решения ОДУ методом Адамса с итерациями
показана на рис. 3.9. Сначала осуществляется ввод (блок 2) зна$
чений первых четырех ординат YÆ, Y1, Y2, Y3 (исходные точки на
рис. 3.9) начального (XÆ) и конечного (XK) значений независи$
Рис. 3.9. Схема алгоритма по методу Адамса
47
мой переменной и шага ее изменений (H), допустимой ошибки EPS
(определяющей конец итерационного процесса), а также макси$
мально допустимого числа итераций KM (обычно KM = 3¸4).
Затем вычисляется многократно используемая в дальнейшем
переменная S (блок 3), значения производных A, B, C, D в соответ$
ствующих начальных точках (блок 4) и исходные значения пере$
менных X и Y (блок 5).
После этого счетчик числа итераций устанавливается в нуль
(блок 6) и вычисляются (блок 7) значения V (части формулы (3.23),
не изменяемой при итерациях) и Y {по формуле (3.22)}.
Затем независимая переменная X получает приращение на за$
данный шаг H (блок 8). Ранее полученное значение Y присваива$
ется переменной W и вычисляется по формуле (3.23) «уточнен$
ное» значение Y (блок 9).
Содержимое счетчика числа итераций увеличивается на едини$
цу (блок 10) и проверяется (блок 11), достигнута ли требуемая
точность при вычислении последних двух итераций.
В случае утвердительного ответа на печать выводятся значения
X, Y, K (блок 12). Если же ответ отрицательный, то перепроверя$
ется (блок 15), не достигнуто ли предельно допустимое число ите$
раций. Если K<KM, то итерационный процесс повторно выполня$
ется переходом к операциям блока 9. Если K = KM (а точность все
еще не достигнута), то следует уменьшить значение шага H. При
этом на печать выводится сообщение «Шаг уменьшить» и выпол$
нение операций алгоритма заканчивается (блок 17).
Если уменьшение шага не требуется, то после вывода на печать
значений X, Y, K (блок 12) проверяется, не достигнута ли верхняя
граница интервала значений независимой переменной X (блок 13).
Если нет (т. е. X < XK), то согласно условиям метода изменяется
значение переменных A, B, C, D (блок 14) с целью продвижения на
один интервал вперед (рис. 3.6). После этого снова выполняются
операции блоков 6–12 и т. д. При X ³ XK вычислительный процесс
завершается.
3.8. Неявные методы численного решения систем ОДУ
При решении на ЭВМ систем ОДУ с плохо обусловленной мат$
рицей [1] часто приходится сталкиваться с необходимостью вос$
производить длительные промежутки времени. В таких случаях
применение перечисленных методов становится малоэффектив$
ным.
48
Подобные явления возникают при численном решении задачи
Коши. Это обстоятельство повлияло на усиленное внимание к разра$
ботке специально машинных методов численного интегрирования,
учитывающих особенности ЭВМ.
Одна из основных идей таких методов заключается в построе$
нии формул численного интегрирования с «обратной связью» – так
называемых неявных методов. В неявных методах можно достичь
увеличения шага h без потери устойчивости вычислительного про$
цесса. В простейшем случае неявный метод получается из форму$
лы (3.19), если воспользоваться формулой правых прямоугольни
ков при взятии интеграла, — yj 11 3 yj 4 hf 1 yj 11 ,tj 11 2.
Аналогично строятся неявные формулы для методов более высо$
кого порядка Эйлера. Например, неявная формула второго порядка
будет иметь вид
1
2 1
2
1
y j 11 5 y j 6 h 3 f tj , y j 6 f tj 11, y j 11 4 ,
8
2 7
а общая формула (3.22) изменится на
1
n
2
y j 11 3 y j 4 h 5 bi y1 j 112i 4 0 hi 11 2 .
i 30
(3.24)
Методы этого класса являются устойчивыми и не имеют ограни$
чений на величину шага дискретности. Шаг дискретности может
выбираться сообразно с требуемой длительностью решения.
Пример 3.3
Применим неявный метод Эйлера для численного интегрирова$
ния уравнения в предыдущем примере. Неявный метод Эйлера име$
ет вид
yk 11 3 yk 4 hf 1 yk 11,t 2.
Применительно к уравнению ( 3.23) при k = n получим:
1
2
yn 21 3 yn 4 h 5104 yn 21 4 e 1hn ,
или
yn 21 3
1
1
3
n
y0
1 4 h10
4
2
n 21
2
4 h 5 1 4 h104 e1
i 1
i
i 112 h
.
(3.25)
49
В этом уравнении условие устойчивости вычислительного процес$
са |1+h•104| >1 не накладывает никаких ограничений на выбор шага
h. Однако точность воспроизведения истинного решения уравнения
(3.23) при большом шаге, h, оказывается низкой, что легко прове$
рить, сопоставив первый член уравнения (3.25) с функцией e–104t, –
истинным решением уравнения (3.23).
При реализации неявных методов на каждом шаге необходимо
решать нелинейные уравнения типа (3.24), что вызывает значитель$
ные осложнения и является недостатком этих методов. Поэтому при$
менение неявных методов сочетают со специальными методами ре$
шения нелинейных уравнений на каждом шаге.
Так, неявный метод Эйлера с модификацией его по формуле Нью$
тона примет вид
11
3
4
5f
yj 21 6 yj 7 91 8 h 1 yj ,tj 21 2 3 8hf 1 yj ,tj 2 4 .
y
5
Напомним, что формула Ньютона для решения уравнения j(y) = 0
имеет вид
11
3 d5 4
yk21 6 yk 7 8
9 5 1 yk 2.
dy yk
Применительно к рассматриваемому случаю в этой формуле сле$
дует положить
yk11 3 yj 11,
yk 3 yj ,
4 1 yk 2 3 y 1 yj 11 2 3 yj 11 5 yj 5 hf 1 yj 11,tj 11 2.
Неявные методы третьей и четвертой степени, получаемые из
(3.22), при соответствующих b i имеют вид
y j 11 3 y j 4
y j 11 3 y j 4
h
15y1 j 11 4 8y1 j 5 y1 j 21 2,
12
h
1 9y1 j 11 4 19y1 j 5 5y1 j 21 4 y1 j 22 2.
24
Численные значения коэффициентов получаются применением
формул, аналогичных (3.20). Из других алгоритмов для неявных
методов приведем методы Куртиса – Хиршфельдера, поскольку их
применение оказалось эффективным при реализации на ЭВМ. Неяв$
ные алгоритмы Куртиса – Хиршфельдера строятся на принципе «диф$
ференцирования назад» и имеют вид:
50
y j 11 1 y j 2 hy1 j 11,
4
1
2
y j 21 3 y j 22 2 hy1 j 11,
3
3
3
8
9
2
6
y j 11 1 y j 21 3 y j 22 2 y j 23 2 hy1 j 11 ,
11
11
11
11
48
36
16
3
12
y j 11 1
y j 21 3
y j 22 2
y j 23 2
y j 24 2
hy1 j 11.
5
25
25
25
25
y j 11 1
(3.26)
Формулы (3.26) строятся на соответствующих аппроксимациях
производной y1 в (3.15)
dy j 11
dt
4
1
y j 11 1 y j 21 dy j 11 y j 11 1 3 y j 21 1 3 y j 22
,
, 1,
2
2
2
h
dt
h
3
т. е. в выполнении своеобразного «дифференцирования назад». Эти
алгоритмы были предложены для интегрирования жестких систем
уравнений Куртисом и Хиршфельдером еще в 1952 г.
Применение неявных методов показывает их преимущества над
явными в смысле устойчивости вычислительного процесса при уве$
личении шага, однако точность получения решений существенно за$
висит от точности решения нелинейных алгебраических уравнений
(3.24). Причем, если в явных методах выбором шага по h удается
получить, хотя и с погрешностью, решение качественно верное, то в
неявных при неправильном выборе шага можно получить устойчи$
вые решения, но качественно отличные от искомого. Значительные
увеличения шага при решении задач неявными методами оказыва$
ются ограниченными из$за резкого снижения точности. Практичес$
ки удается увеличить шаг не более чем в 5–10 раз.
Разностное уравнение, заменяющее исходное (3.4), является обоб$
щенной формулой методов численного интегрирования, распростра$
ненных в вычислительной практике:
r
9 37ak zn1k 5 hbkf 1 tn1k ,zm1k 248 6 0.
n 20
(3.27)
Из этой формулы в зависимости от вида f и значений параметров
ak, bk, r получаются как классические методы (3.7), (3.24), так и
разработанные специально для ЭВМ при современной их классифи$
кации. Так, если r = 1, то из (3.27) получаются одношаговые методы
51
(см., например, формулы Рунге – Кутты), если r > 1 – многошаговые
{cм. (3.24)}, если br = 0, то получаем явные методы (3.22); если b2 ¹ 0
– неявные (3.24). В последнее время разработан ряд вычислитель$
ных процедур на основе явных и неявных методов численного интег$
рирования.
Эффективное применение неявные методы и алгоритмы нашли в
процедуре, предложенной Гиром. В процедуре Гира применяются ме$
тоды (3.24), (3.26). Однако это не исключает возможности примене$
ния и других численных методов. Неявный характер употребляемых
в процедуре Гира методов предполагает решение на каждом шаге си$
стемы алгебраических уравнений, для чего предусматривается ис$
пользование трех модификаций итерационного метода Ньютона.
Для этого алгоритмы (3.24) представляются в виде
k1
k2
j 20
j 20
yn 1 5 2 j yn 1 j 3 h 5 4j y1 n 1 j ,
(3.28)
где aj и bj – числовые коэффициенты, а k1 k2 – количество соответству$
ющих слагаемых в (3.28), a = 1.
Решение нелинейных алгебраических уравнений относительно
(3.28) строится как следующий процесс:
1k2
1 k 112
yn 5 yn
1 k2
6 P 21 * 7 38 yn 49,
k
k
1
2
1k2
1 k2
1k2
1k2
где 5 3
yn 4 6 yn 7 h80 y1 n 7 9 j y1 n 1 j 7 h 8i y1 n 1 j ,
j 21
j 21
P5
34
3 yn
5 1 6 h70
1
k
3f yn1 2 ,tn
1k2
yn 1 y n
k
3yn1 2
2,
где ynk – kе приближение yk в итерационном процессе k Î[1, N], а f(yn,
t n) = y'n.
В случае системы уравнений (3.3), вместо скалярной величины p
строится матрица
P4
3F
3 Yn
4 E 5 h60
1k2
Y n 1Y n
1
3F t0 , Ynk
3Y
1k2
2.
В процедуре Гира автоматически выбирается численный метод и
шаг по заданной точности вычислений.
52
Если в некоторый момент времени используется метод степени l
(например, в (3.28) k1 = l, k2 = l), то одновременно оценивается ло$
кальная погрешность методов степени 1 и 2. В дальнейшем применя$
ется тот метод, который при заданной точности позволяет интегри$
ровать с большим значением шага h.
Еще одно преимущество процедуры Гира в том, что она не требует
стартового алгоритма. Интегрирование начинается с методов первой
степени, а затем степень метода при необходимости изменяется в со$
ответствии с задаваемой пользователем точностью. Следует отме$
тить, что алгоритмы процедуры подвержены неустойчивости, кото$
рая проявляется в колебательном характере шага, принимающем то
разумные значения, то неприемлемо малые. Возникающие трудно$
сти можно исключить, уточняя соответствующие константы aj и bj
эмпирическим путем.
Все рассмотренные методы являются скалярными, т. е. для каж$
дого уравнения системы (3.3) составляется соответствующее разно$
стное уравнение. Поэтому они обладают присущими скалярным ме$
тодам недостатками при применении их к системам большой размер$
ности. Качественно новым является создание системных (матрич$
ных) методов, которые позволяют учитывать свойства правой части
системы (3.3) и обеспечивают независимость выбора шага от устой$
чивости вычислительного процесса.
3.9. Задания для выполнения лабораторной работы № 3
Найти численное решение y(x) обыкновенного дифференциально$
го уравнения при заданных начальных условиях на отрезке [a, b] с
шагом h, используя заданный метод численного интегрирования
(табл. 3.1).
Найти аналитическое решение (вручную). В программе предус$
мотреть вычисление отклонения d численного решения от анали$
тического в точках вывода (шаг вывода взять равным h). Исследо$
вать связь между указанным отклонением и шагом интегрирова$
ния, для чего провести вычисления повторно при шаге h/4 и срав$
нить полученное отклонение *d с отклонением d. Результаты печа$
тать в виде таблиц значений x, y, d при шаге h и при шаге h/4 [7].
Эти методы будут изложены во второй части методических указа$
ний.
53
Таблица 3.3
54
Библиографический список
1. Сольницев Р. И. Автоматизация проектирования систем управле$
ния, гл. 5. М.: ВШ, 1991
2. Сольницев Р. И. , Гришанова Л. И., Слюсаренко А. С. Математическое
обеспечение информационных технологий. Часть I. Непрерывные систе$
мы. СПб, ГУАП, 2004.
3. Турчак Л. И. Основы численных методов. М.: Наука, 1987.
4. Демидович Б. П., Марон И. А. Основы вычислительной математики.
М.: Наука, 1966.
5. Алексеев А. В., Ковтун И. В., Потоцкий В.К., Пресняк А. С., Тертеро
ва И. М., Тревогина И. В., Щадилов А. Е. Вычислительная математика и
программирование: Метод. указания к выполнению лабораторных работ/
ЛИАП. Л., 1988.
Оглавление
Предисловие .......................................................................
3
Лабораторная работа № 1 ....................................................
4
ЧИСЛЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНЫХ УРАВНЕНИЙ .......... 4
1.1. Общие методические указания к численным методам решения
нелинейных уравнений .....................................................
4
1.2. Метод половинного деления .........................................
6
1.3. Метод последовательных приближений .........................
8
1.4. Метод последовательных приближений с переходом
к обратным функциям ......................................................
11
1.5. Метод касательных (метод Ньютона) .............................
12
1.6. Упрощенный метод касательных ..................................
14
1.7. Метод хорд ................................................................
15
1.8. Модифицированный метод хорд ...................................
18
1.9. Комбинированный метод (вариант 1) .............................
19
1.10. Комбинированный метод (вариант 2) ...........................
20
1.11. Комбинированный метод (вариант 3) ...........................
22
1.12. Метод последовательных приближений
для системы двух нелинейных уравнений ...........................
23
1.13. Методические указания к составлению алгоритма решения
нелинейного уравнения (системы уравнений) .......................
24
1.14. Варианты заданий для выполнения лабораторной работы ....
25
Лабораторная работа№ 2 ..................................................... 2 7
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ........................................ 2 7
2.1. Постановка задачи ......................................................
2.2. Метод прямоугольников ..............................................
2.3. Формула трапеций .....................................................
2.4. Метод Симпсона .........................................................
2.5. Задания для выполнения лабораторной работы ...............
27
27
28
29
29
Лабораторная работа № 3 .................................................... 3 1
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ........................... 3 1
3.1. Общие сведения о дифференциальных уравнениях .........
3.2. Классификация численных методов решения
обыкновенных дифференциальных уравнений .....................
3.3. Метод степенного ряда ................................................
3.4. Метод Эйлера и его модификации ................................
3.5. Метод Рунге – Кутты ..................................................
3.6. Метод Адамса ............................................................
3.7. Алгоритмизация методов численного решения ОДУ .......
3.8. Неявные методы численного решения систем ОДУ .........
3.9. Задания для выполнения лабораторной работы № 3 .........
Библиографический список ...............................................
31
32
40
41
43
44
45
48
53
55
55
Документ
Категория
Без категории
Просмотров
0
Размер файла
728 Кб
Теги
grishanovasokmet
1/--страниц
Пожаловаться на содержимое документа