Рынок: модели, черты
Автор работы: Пользователь скрыл имя, 01 Декабря 2013 в 18:54, реферат
Краткое описание
Инженерные и научные задачи часто приводят к решению
различных уравнений или систем уравнений, описывающих по-
ведение параметров объекта, например динамические нагрузки
на строительную конструкцию или тепловые потоки через сте-
ны дома. Совокупность всех уравнений и дополнительных усло-
вий, которым должно удовлетворять решение, называется мате-
матической моделью.
Содержание
ВВЕДЕНИЕ.................................................................................................. 4
1. ПРИБЛИЖЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ........................................................ 8
1.1. П
ОСТАНОВКА ЗАДАЧИ
....................................................................... 8
1.2. П
РИБЛИЖЕННЫЕ МЕТОДЫ
.................................................................. 9
1.3. С
ТАНДАРТНЫЕ ФУНКЦИИ
M
ATH
CAD............................................ 16
2. МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ
УРАВНЕНИЙ.............................................................................................21
2.1. О
БЩИЕ ВОПРОСЫ
..............................................................................21
2.2. Т
ОЧНЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
............................................................... 22
2.3. И
ТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ ЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ СИСТЕМ
..................................................................... 27
2.4. С
ТАНДАРТНЫЕ ФУНКЦИИ ПАКЕТА
M
ATH
CAD................................ 31
3. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ......................... 33
3.1. П
ОСТАНОВКА ЗАДАЧИ ИНТЕРПОЛЯЦИИ
........................................... 33
3.2. Л
ОКАЛЬНАЯ ИНТЕРПОЛЯЦИЯ
........................................................... 34
3.3. Г
ЛОБАЛЬНАЯ ИНТЕРПОЛЯЦИЯ
.......................................................... 38
3.4. П
ОЛИНОМ
Л
АГРАНЖА
...................................................................... 39
3.5. М
ЕТОД НАИМЕНЬШИХ КВАДРАТОВ
................................................. 41
4. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ... 45
4.1. Ч
ИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ
................................................ 45
4.2. И
СПОЛЬЗОВАНИЕ СТАНДАРТНЫХ ФУНКЦИЙ
M
ATH
CAD
ДЛЯ ДИФФЕРЕНЦИРОВАНИЯ
................................................................... 48
4.3. Ч
ИСЛЕННОЕ ИНТЕГРИРОВАНИЕ
....................................................... 51
4.4. И
СПОЛЬЗОВАНИЕ СТАНДАРТНЫХ ФУНКЦИЙ
M
ATH
CAD
ДЛЯ ИНТЕГРИРОВАНИЯ
............................................................................ 54
5. ЧИСЛЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ.............................................. 59
5.1. Ч
ИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ
К
ОШИ
............................. 59
5.2. К
РАЕВАЯ ЗАДАЧА ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ
УРАВНЕНИЙ ВТОРОГО ПОРЯДКА
.............................................................. 63
5.3. Ж
ЕСТКИЕ ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ
...... 66
5.4. Р
ЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
И СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
В ПАКЕТЕ
M
ATH
CAD............................................................................ 69
6. РЕШЕНИЕ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ........... 75
6.1. О
СНОВНЫЕ ПОНЯТИЯ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ
..... 75
6.2. П
АРАБОЛИЧЕСКИЕ УРАВНЕНИЯ
....................................................... 77
6.3. П
РИБЛИЖЕННЫЕ МЕТОДЫ РЕШЕНИЯ УРАВНЕНИЙ
ГИПЕРБОЛИЧЕСКОГО ТИПА
..................................................................... 89
6.4. П
РИБЛИЖЕННЫЕ МЕТОДЫ РЕШЕНИЯ УРАВНЕНИЯ
П
УАССОНА
....... 94
ЗАКЛЮЧЕНИЕ......................................................................................... 97
БИБЛИОГРАФИЧЕСКИЙ СПИСОК..................................................... 98
Прикрепленные файлы: 1 файл
3
ЧИСЛЕННЫЕ МЕТОДЫ
РЕШЕНИЯ ИНЖЕНЕРНЫХ ЗАДАЧ
В ПАКЕТЕ MathCAD
y
y=F(x)
x
x
0
x
1
x
2
x
3
x*
НОВОСИБИРСК 2005
f x()
x
2
x
−
1
−
:=
df x()
f x()
d
dx
:=
newt x ε,
( )
x
x
f x
( )
df x
( )
−
←
f x()
ε
>
while
x
:=
newt 0 0.01
,
(
)
0.619
−
=
И.А. БЕДАРЕВ
О.Н. БЕЛОУСОВА
Н.Н. ФЕДОРОВА
4
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ
РОССИЙСКОЙ ФЕДЕРАЦИИ
НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ
АРХИТЕКТУРНО-СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ (СИБСТРИН)
И.А. Бедарев, О.Н. Белоусова, Н.Н. Федорова
ЧИСЛЕННЫЕ МЕТОДЫ
РЕШЕНИЯ ИНЖЕНЕРНЫХ ЗАДАЧ
В ПАКЕТЕ MathCAD
УЧЕБНОЕ ПОСОБИЕ
НОВОСИБИРСК 2005
5
ОГЛАВЛЕНИЕ
ВВЕДЕНИЕ.................................................................................................. 4
1. ПРИБЛИЖЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ........................................................ 8
1.1. П
ОСТАНОВКА ЗАДАЧИ
....................................................................... 8
1.2. П
РИБЛИЖЕННЫЕ МЕТОДЫ
.................................................................. 9
1.3. С
ТАНДАРТНЫЕ ФУНКЦИИ
M
ATH
CAD............................................ 16
2. МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ
УРАВНЕНИЙ.............................................................................................21
2.1. О
БЩИЕ ВОПРОСЫ
..............................................................................21
2.2. Т
ОЧНЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
............................................................... 22
2.3. И
ТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ ЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ СИСТЕМ
..................................................................... 27
2.4. С
ТАНДАРТНЫЕ ФУНКЦИИ ПАКЕТА
M
ATH
CAD................................ 31
3. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ......................... 33
3.1. П
ОСТАНОВКА ЗАДАЧИ ИНТЕРПОЛЯЦИИ
........................................... 33
3.2. Л
ОКАЛЬНАЯ ИНТЕРПОЛЯЦИЯ
........................................................... 34
3.3. Г
ЛОБАЛЬНАЯ ИНТЕРПОЛЯЦИЯ
.......................................................... 38
3.4. П
ОЛИНОМ
Л
АГРАНЖА
...................................................................... 39
3.5. М
ЕТОД НАИМЕНЬШИХ КВАДРАТОВ
................................................. 41
4. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ... 45
4.1. Ч
ИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ
................................................ 45
4.2. И
СПОЛЬЗОВАНИЕ СТАНДАРТНЫХ ФУНКЦИЙ
M
ATH
CAD
ДЛЯ ДИФФЕРЕНЦИРОВАНИЯ
................................................................... 48
4.3. Ч
ИСЛЕННОЕ ИНТЕГРИРОВАНИЕ
....................................................... 51
4.4. И
СПОЛЬЗОВАНИЕ СТАНДАРТНЫХ ФУНКЦИЙ
M
ATH
CAD
ДЛЯ ИНТЕГРИРОВАНИЯ
............................................................................ 54
5. ЧИСЛЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ.............................................. 59
5.1. Ч
ИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ
К
ОШИ
............................. 59
5.2. К
РАЕВАЯ ЗАДАЧА ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ
УРАВНЕНИЙ ВТОРОГО ПОРЯДКА
.............................................................. 63
5.3. Ж
ЕСТКИЕ ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ
...... 66
5.4. Р
ЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
И СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
В ПАКЕТЕ
M
ATH
CAD............................................................................ 69
6. РЕШЕНИЕ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ........... 75
6.1. О
СНОВНЫЕ ПОНЯТИЯ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ
..... 75
6.2. П
АРАБОЛИЧЕСКИЕ УРАВНЕНИЯ
....................................................... 77
6.3. П
РИБЛИЖЕННЫЕ МЕТОДЫ РЕШЕНИЯ УРАВНЕНИЙ
ГИПЕРБОЛИЧЕСКОГО ТИПА
..................................................................... 89
6.4. П
РИБЛИЖЕННЫЕ МЕТОДЫ РЕШЕНИЯ УРАВНЕНИЯ
П
УАССОНА
....... 94
ЗАКЛЮЧЕНИЕ......................................................................................... 97
БИБЛИОГРАФИЧЕСКИЙ СПИСОК..................................................... 98
6
ВВЕДЕНИЕ
Инженерные и научные задачи часто приводят к решению
различных уравнений или систем уравнений, описывающих по-
ведение параметров объекта, например динамические нагрузки
на строительную конструкцию или тепловые потоки через сте-
ны дома. Совокупность всех уравнений и дополнительных усло-
вий, которым должно удовлетворять решение, называется мате-
матической моделью. Простая математическая модель – это со-
вокупность алгебраических формул, по которым явно вычисля-
ются искомые величины. Однако чаще всего поведение пара-
метров описывается дифференциальными уравнениями в част-
ных производных. Найти решение этих сложных задач можно
только с использованием современных быстродействующих
ЭВМ. Решение сложной математической задачи на ЭВМ вклю-
чает в себя необходимые этапы выбора метода решения, созда-
ния алгоритма, разработки программы и ее тестирования. После
этого можно применять разработанный пакет программ для ре-
шения нужной задачи. Даже для того, чтобы воспользоваться
стандартной, т.е. уже готовой программой, нужно иметь пред-
ставление о существующих методах решения, их преимущест-
вах, недостатках и особенностях использования.
Все математические задачи классифицированы, т.е. объе-
динены в некоторые группы. Для каждой группы задач сущест-
вует набор стандартных методов, которые изучает специальный
раздел математики – «Вычислительная математика» или «Мето-
ды вычислений».
Все методы решения уравнений можно разделить на два
класса: точные и приближенные. В точных методах решение по-
лучают в виде формул за конечное число операций, но их можно
применять только для решения уравнений специального вида.
В общем случае задачу можно решить только приближенно.
Приближенные методы позволяют получить решение в виде
бесконечной последовательности, сходящейся к точному ре-
шению.
Использование ЭВМ выдвигает дополнительные требова-
ния к алгоритму нахождения как точного, так и приближенного
7
решения: он должен быть устойчивым, реализуемым и эконо-
мичным. Устойчивость означает, что малые погрешности, вне-
сенные в процессе решения, не приводят к большим ошибкам в
конечном результате. Погрешности возникают из-за неточного
задания исходных данных (неустранимые ошибки), из-за округ-
ления чисел, которое всегда имеет место при расчетах на ЭВМ,
а также связаны с точностью используемого метода. Реализуе-
мость алгоритма означает, что решение может быть получено
за допустимое время. При этом надо иметь в виду, что время
приближенного решения зависит от точности, с которой мы хо-
тим получить решение. На практике точность выбирают с уче-
том реализуемости алгоритма на той ЭВМ, которую предпола-
гается использовать для решения. Экономичным называется ал-
горитм, который позволяет получить решение с заданной точно-
стью за минимальное количество операций и, следовательно, за
минимальное расчетное время.
В изучаемом курсе мы познакомимся с основными метода-
ми, применяемыми для решения различных математических за-
дач. Первым рассматриваемым классом задач будут нелинейные
алгебраические уравнения. Потом мы научимся решать системы
линейных алгебраических уравнений и обыкновенные диффе-
ренциальные уравнения, приближенно находить производные и
интегралы, а также познакомимся с основными понятиями ин-
терполяции (приближения) функций. Заключительная глава по-
священа приближенному решению уравнений в частных произ-
водных.
Каждая тема, кроме теоретического материала, содержит
примеры использования методов для решения конкретных за-
дач, описания основных вычислительных алгоритмов, тексты
программ и описание стандартных функций пакета MathCAD,
реализующих изученные вычислительные алгоритмы.
8
1. ПРИБЛИЖЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
1.1. Постановка задачи
Дано нелинейное алгебраическое уравнение
F(x) = 0.
(1.1)
Нелинейность уравнения означает, что график функции не есть
прямая линия. Решить уравнение – это найти x* ⊂ R : F(x*) = 0.
y
y =F(x)
x
x
1
*
x
2
*
x
3
*
Значение x* называют кор-
нем уравнения. Нелинейное
уравнение может иметь не-
сколько корней. Геометри-
ческая интерпретация такой
ситуации представлена на
рис. 1.1. Корнями уравнения
(1.1) являются точки x
1
*, x
2
*,
x
3
*, в которых функция F(x)
пересекает ось x.
Рис. 1.1. Геометрическая иллюстрация
уравнения (1.1)
Необходимое условие существования корня уравнения (1.1)
и достаточное условие единственности следуют из известной
теоремы Больцано–Коши. Пусть F(x) непрерывна и F(a)F(b) < 0
(т.е. на концах интервала функция имеет разные знаки). Тогда
внутри отрезка [a, b] существует корень уравнения F(x) = 0. Ко-
рень будет единственным, если F
′
(x) не меняет знак на отрезке
[a, b], т.е. F(x) – монотонная функция.
Методы решения уравнения (1.1) можно разделить на точ-
ные (аналитические) и приближенные (итерационные). Точными
методами корень находится за конечное число действий и пред-
ставляется некоторой алгебраической формулой. Процесс нахо-
ждения решения приближенными методами бесконечен. Реше-
нием называется бесконечная последовательность {x
n
}, такая,
что
x*
x
n
n
=
∞
→
lim
. По определению предела, для любого сколь
угодно малого наперед заданного
ε
найдется такое N, что при
n > N |x
n
– x*| <
ε
. Члены этой последовательности {x
n
} называ-
ются последовательными приближениями к решению, или ите-
рациями. Наперед заданное число
ε
называют точностью мето-
9
да, а N – это количество итераций, которое необходимо выпол-
нить, чтобы получить решение с точностью
ε
. Существуют раз-
личные методы нахождения приближенного решения, т.е. спо-
собы построения последовательности итераций {x
n
}, однако все
они имеют общие этапы, изображенные на рис. 1.2.
Используются различ-
ные критерии остановки
итерационного процесса:
– |x
n
– x*| <
ε
. К сожале-
нию, это условие не все-
гда возможно проверить,
т.к. x* неизвестно;
– ⎟F(x
n
)⎟ <
ε
, где F(x
n
) – не-
вязка метода;
– |x
n+1
– x
n
| <
ε
, т.е. разница
между соседними итера-
циями стала мала.
Задать начальное
приближение x
0
Проверить
условие выхода из
итерационного
процесса
Найти следующее
приближение
x
n+1
=
ϕ
(x
n
, x
n–1
, …, x
1
, x
0
)
нет
да
Рис. 1.2. Этапы итерационного
процесса
1.2. Приближенные методы
Прежде чем использовать приближенный метод, надо ис-
следовать уравнение на наличие корней и уточнить, где эти кор-
ни находятся, т.е. найти интервалы изоляции корней. Интерва-
лом изоляции корня называется отрезок, на котором корень
уравнения существует и единственный. Каждому корню соот-
ветствует свой интервал изоляции. Если корней несколько, то
для каждого нужно найти интервал изоляции. Существуют раз-
личные способы исследования функции: аналитический, таб-
личный, графический. Аналитический способ состоит в нахож-
дении экстремумов функции F(x), исследовании ее поведения
при
,∞
→
x
нахождении участков возрастания и убывания
функции. Табличный способ – это построение таблицы, состоя-
щей из столбца аргумента x и столбца значений функции F(x).
О наличии корней свидетельствуют перемены знака функции.
Чтобы не произошло потери корней, шаг изменения аргумента
должен быть достаточно мелким, а интервал изменения
–
доста-
10
точно большим. Графический способ – это построение графика
функции F(x) и определение числа корней по количеству пере-
сечений графика с осью x. Ниже для иллюстрации графическим
способом исследовано уравнение F(x) = 0.4 ⋅ 2
x
– 0.5x – 1 = 0.
4
2
0
2
0
2
4
3.2
1
−
f1 x
( )
f2 x
( )
3
5
−
x
На рис. 1.3 приведены по-
строенные
с
помощью
MathCAD графики функций
F
1
(x) = 0.4 ⋅ 2
x
и F
2
(x) = 0.5x
+ 1. Корнями являются точ-
ки, в которых пересекаются
два графика. Рисунок пока-
зывает, что исходное урав-
нение имеет два корня, рас-
положенных на интервалах
[
–
3, 0] и [0, 3].
Рис. 1.3. Графический способ нахо-
ждения интервалов изоляции
Пусть интервалы изоляции корней известны. Познакомимся
с несколькими итерационными методами, позволяющими найти
корень на известном интервале изоляции [a, b].
Метод деления отрезка пополам (дихотомии). Найдем сере-
дину отрезка [a, b]: c = (a+b)/2. Корень остался на одной из час-
тей: [a, c] или [c, b]. Если F(a)
⋅
F(с) < 0, то корень попал на отре-
зок [a, c], тогда деление отрезка можно повторить, приняв в каче-
стве нового правого конца точку c, т.е. b = c. В противном случае
корень попал на половину [c, b], и необходимо изменить значение
левого конца отрезка: a = c. Поскольку корень всегда заключен
внутри отрезка, итерационный процесс можно останавливать, ес-
ли длина отрезка станет меньше заданной точности: |b – a| <
ε
.
В случае сложных уравнений вычисления необходимо про-
водить с использованием ЭВМ. На рис. 1.4 приведен текст про-
граммы MathCAD, реализующей метод дихотомии для решения
уравнения F(x) = x
2
–
x – 1 = 0. Метод реализован в виде функ-
ции от аргументов a, b (концы интервала изоляции) и точности
метода
ε
. Вызов этой функции при значениях a = –1, b = 0,
ε
= 0.01 дает значение корня – 0.617, при этом невязка уравне-
ния, которую вычисляют, чтобы убедиться в правильности ре-
шения, равна –1.892⋅10
–3
.
11
f x()
x
2
x
−
1
−
:=
mdp a b, ε,
(
)
c
a
b
+
(
)
2
←
a
c
←
f a() f c()
⋅
0
>
if
b
c
←
otherwise
b a
−
ε
>
while
c
:=
mdp 1
− 0, 0.01
,
(
) =
f mdp 1
− 0, 0.01
,
(
)
(
) =
Рис. 1.4. Метод деления отрезка пополам
Метод хорд. В этом методе кривая F(x) заменяется прямой
линией – хордой, стягивающей точки (a, F(a)) и (b, F(b)). В за-
висимости от знака выражения F(a)F
″
(a) метод хорд имеет два
варианта, изображенных на рис. 1.5.
Пусть F(a)F
″
(a)>0 (см. рис. 1.5а). Тогда x
0
= b, точка a будет
оставаться неподвижной. Следующее приближение x
1
находим как
точку пересечения хорды, соединяющей точки (a, F(a)) и (x
0
, F(x
0
))
с осью x. Уравнение хорды:
)
(
)
(
)
(
)
(
0
0
a
x
a
x
a
F
x
F
a
F
y
−
−
−
+
=
.
Тогда точка пересечения хорды с осью x:
)
(
)
(
)
)(
(
0
0
1
a
F
x
F
a
x
a
F
a
x
−
−
−
=
.
Пусть теперь F(a)F
″
(a) < 0 (рис. 1.5б). Тогда x
0
= a, точка b
неподвижна. Проведем хорду, соединяющую точки (b, F(b)) и
(x
0
, F(x
0
)):
)
(
)
(
)
)
(
0
0
0
0
x
x
x
b
x
F
F(b
x
F
y
−
−
−
+
=
. Вычисляем точку
пересечения хорды с осью x:
)
(
)
(
)
)(
(
0
0
0
0
1
x
F
b
F
x
b
x
F
x
x
−
−
−
=
.
На следующей итерации в качестве x
0
надо взять вычислен-
ное значение x
1
. Окончание итерационного цикла в этом методе
происходит по условию малости невязки уравнения: |F(x
1
)| <
ε
.
12
y
y = F(x)
x
x*
x
1
x
2
a
x
0
= b
а
y
y = F(x)
x
x*
x
1
x
2
x
0
=a
b
б
Рис. 1.5. Метод хорд: а – F(a)F″(a) > 0; б – F(a)F″(a) < 0
Метод Ньютона (касательных). Как и предыдущий, этот
метод основан на замене исходного нелинейного уравнения (1.1)
линейным уравнением, которое можно легко решить. Иллюст-
рация метода представлена на рис. 1.6. Пусть x
0
– начальное
приближение. Построим касательную к функции y = F(x), про-
ходящую через точку (x
0
, F(x
0
)). Найдем пересечение касатель-
ной:
)
)(
(
)
(
0
0
0
x
x
x
F
x
F
y
−
′
+
=
y
y =F(x)
x
x
0
x
1
x
2
x
3
x*
с осью x:
)
(
)
(
0
0
0
1
x
F
x
F
x
x
′
−
=
.
На следующей итерации в
качестве x
0
надо взять вы-
численное
значение
x
1
.
Окончание итерационного
цикла, как и в методе хорд,
выполняется по невязке
уравнения: |F(x
1
)| <
ε
.
Рис. 1.6. Метод Ньютона
Как показывают практика и теоретические оценки, метод
Ньютона позволяет достаточно быстро получить решение. Не-
достатком метода является то, что для некоторых функций F(x)
при неудачном выборе начального приближения метод расходит-
ся. Ситуацию легко исправить, если выбрать x
0
ближе к x*.
Упрощенный метод Ньютона. Эта модификация метода
Ньютона используется, если производная F
′
(x) представляет со-
бой сложную функцию и для ее вычисления на каждой итерации
13
требуется много времени. Зададим x
0
– начальное приближение
и вычислим производную z = F
′
(x
0
). На следующих итерациях
используется
вычисленное
значение
производной:
z
x
F
x
x
k
k
k
)
(
1
−
=
+
. Это упрощение несколько замедляет процесс
сходимости к решению, однако сокращает время каждого ите-
рационного цикла.
Метод Чебышева является развитием метода Ньютона. Если
приближение x
k
известно, то следующее приближение находим
по формуле
3
1
))
(
(2
)
(
)
(
)
(
)
(
k
k
k
k
k
k
k
x
f
x
f
x
f
x
f
x
f
x
x
′
′′
−
′
−
=
+
. Этот метод позво-
ляет быстро получить корень, однако требует еще более точного
задания начального приближения x
0
.
Метод секущих применяется, когда вычисление производной
F
′
(x) занимает много времени, а также если функция F(x) задана
таблично. Для этого метода необходимо задать два начальных
приближения: x
0
, x
1
, поэтому он называется двухшаговым (все рас-
смотренные выше методы были одношаговыми). Значение произ-
водной можно вычислить с помощью конечно-разностного соот-
ношения F
′
(x
1
)
≈
0
1
0
1
)
(
)
(
x
x
x
F
x
F
−
−
. Подставляя это соотношение в
формулу для метода Ньютона, получим
)
(
)
(
)
)(
(
0
1
0
1
1
1
2
x
F
x
F
x
x
x
F
x
x
−
−
−
=
.
На следующей итерации в качестве приближений x
0
, x
1
возьмем
уже вычисленные значения x
1
, x
2
, т.е. x
0
= x
1
, x
1
= x
2
, и будем вычис-
лять новое приближение x
2
. Окончание итерационного цикла про-
изводится по невязке уравнения: |F(x
2
)| <
ε.
Метод простой итерации (МПИ). Этот метод является
обобщением всех описанных выше одношаговых методов. Слово
«простой» означает, что для вычисления следующего приближе-
ния необходимо знать только одно предыдущее приближение. С
помощью эквивалентных преобразований приведем исходное
уравнение (1.1) к виду, удобному для применения метода простой
итерации: x =
ϕ
(x). Выберем начальное приближение x
0
∈ [a, b].
14
Следующие итерации находим по формуле: x
k+1
=
ϕ
(x
k
). Итера-
ционный процесс заканчивается, если |
ϕ
(x
k
) – x
k
| <
ε
или, что
то же самое, |x
k+1
– x
k
| <
ε
. Иллюстрация метода представлена на
рис. 1.7.
y
y = x
x
x*
x
0
x
1
ϕ
(x
0
)
y =
ϕ
(x)
x
2
ϕ
(x
1
)
а
y = x
x
x*
x
0
x
1
y =
ϕ
(x)
x
3
x
2
x
4
б
y
Рис. 1.7. Сходящийся метод простой итерации:
а –
ϕ′
(x) > 0;
б –
ϕ′
(x) < 0
Способ сведения исходного уравнения к виду x =
ϕ
(x) очень
важен, поскольку от вида функции
ϕ
(x) зависит, будет ли итера-
ционный процесс сходиться или нет.
Чтобы понять, каким условиям должна удовлетворять поро-
ждающая итерации функция
ϕ
(x), проведем некоторые теорети-
ческие оценки. Пусть x
k
– известное приближение, отличающее-
ся от искомого корня на величину погрешности |x
k
–
x*|. Сле-
дующее приближение x
k+1
будет отличаться от корня на величи-
ну |x
k+1
– x*| = |
ϕ
(x
k
) – x*| = |
ϕ
(x
k
) –
ϕ
(x*)| = |
ϕ′
(
ξ
)| |x
k
– x*|, где
ξ
–
некоторая средняя точка отрезка (x
k
, x*). Здесь использована
теорема Лагранжа о конечном приращении, а также равенство
x* =
ϕ
(x*), которое справедливо, поскольку x* – корень. В схо-
дящихся методах каждое следующее приближение должно быть
ближе к корню, чем предыдущее, а это справедливо, если
|
ϕ′
(
ξ
)| ≤ q < 1. Понятно, что чем меньше число q, тем быстрее
сходится итерационный процесс. Можно получить оценку:
|x
k+1
– x*| < q⋅|x
k
– x*| < q
2
⋅|x
k–1
– x*| < … < q
k
⋅|x
0
– x*|.
Если удастся показать, что |x
k+1
– x*| < q⋅|x
k
– x*|
α
,
α
> 1, то гово-
рят, что метод сходится с порядком
α
.
Можно доказать, что необходимым условием сходимости
МПИ является |
ϕ′
(x*)| < 1. Это условие трудно проверить, т.к.
корень неизвестен. Но можно проверить более сильное доста-
15
точное условие сходимости МПИ |
ϕ′
(x)| < 1, где x – любая точка
некоторого отрезка [a, b], содержащего корень.
Проверим, выполняется ли необходимое условие метода
Ньютона, для которого
)
(
)
(
)
(
x
F
x
F
x
x
′
−
=
ϕ
. Вычислим производ-
ную:
(
)
(
)
(
)
1
0
)*
(
)*
(
)*
(
)*
(
)*
(
)*
(
)*
(
1
)*
(
2
2
2
<
=
′
′′
=
′
′′
−
′
−
=
′
x
F
x
F
x
F
x
F
x
F
x
F
x
F
x
ϕ
.
Следовательно, метод Ньютона всегда сходится в некоторой ок-
рестности корня x*. Можно показать, что |x
k+1
– x*|<q⋅|x
k
– x*|
2
,
т.е.
α
= 2: метод сходится со вторым порядком. Порядок сходи-
мости метода Чебышева еще выше (
α
= 3), для метода секущих
α
= 1.67, а для метода дихотомии
α
= 1.
Если необходимое условие нарушается, то МПИ будет рас-
ходиться. Иллюстрация расходящихся итерационных методов
приведена на рис. 1.8.
y
y = x
x
x
0
x
1
y = ϕ(x)
x
2
а
x*
y=x
x
x*
x
0
x
1
y =
ϕ
(x)
x
3
x
2
x
4
б
y
Рис. 1.8. Расходящийся метод простой итерации:
а –
ϕ′
(x) > 0;
б –
ϕ′
(x) < 0
Метод релаксации – это универсальный вариант МПИ, в ко-
тором
ϕ
(x) = x –
τ
F(x). Параметр релаксации
τ ≠
0 подберем таким
образом, чтобы выполнялось достаточное условие сходимости
метода: ⎟
ϕ′
(x)⎟ < 1 для всех x ∈ [a, b]. В методе релаксации
⎟
ϕ′
(x)⎟ = |1 –
τ⋅
F
′
(x)| < 1, тогда –2 <
τ⋅
F
′
(x) < 0. Таким образом, при
F
′
(x)>0, x ∈ [a, b], условием сходимости будет –2/F
′
(x)<
τ
< 0.
При F
′
(x) < 0 параметр надо выбирать из условия –2/F
′
(x)>
τ
> 0.
16
Ясно, что при малых
τ
МПИ будет сходиться медленно, т.к. рас-
стояние между соседними итерациями |x
k+1
– x
k
| = |
τ
|
⋅
|F
′
(x)| ма-
ло. Можно ли выбрать такое значение параметра, при котором
скорость сходимости будет максимальной? Несложный анализ
показывает, что оптимальным значением будет
τ
= 2/(M + m),
где
),
(
max
]
,
[
x
F
M
b
a
x
′
=
∈
),
(
min
]
,
[
x
F
m
b
a
x
′
=
∈
а скорость сходимости в
этом случае будет определяться константой
m
M
m
M
q
+
−
=
.
1.3. Стандартные функции MathCAD
Для решения уравнения (1.1) в MathCAD служит функция
root, реализующая описанный выше метод секущих. Если
F(x) – это полином, то вычислить все его корни можно также с
помощью функции polyroots.
Встроенная функция root в зависимости от типа задачи мо-
жет иметь либо два аргумента: root (f(x), x), либо четыре аргу-
мента: root (f(x), x, a, b). Здесь f(x) – скалярная функция, опре-
деляющая уравнение (1.1); x – скалярная переменная, относи-
тельно которой решается уравнение; a, b – границы интервала,
внутри которого происходит поиск корня.
Первый тип функции root требует дополнительного зада-
ния начального значения переменной x. Для этого нужно просто
предварительно присвоить x некоторое число. Поиск корня бу-
дет производиться вблизи этого числа. Таким образом, присвое-
ние начального значения требует априорной информации о
примерной локализации корня.
Рассмотрим решение уравнения sin(x) = 0, которое имеет
бесконечное количество корней x
N
= N
π
(N = 0, ±1, ±2, ...). Для
поиска корня средствами MathCAD требуется его предвари-
тельная локализация путем задания начального приближения,
например, x = 0.5. MathCAD находит с заданной точностью
только один корень x
0
= 0, лежащий наиболее близко к заданно-
му начальному приближению. Если задать другое начальное
значение, например, x = 3, то решением будет другой корень
уравнения x
1
= π и т.д.
17
На рис. 1.9 приведен пример вызова стандартной функции
root с двумя аргументами для нахождения корней уравнения
sin(x) = 0, график функции f(x) = sin(x) и положение найденного
корня.
x 0.5
:=
f x()
sin x()
:=
s
root f x() x,
(
)
:=
s
6.2
−
10
7
−
×
=
1
0
1
1
0
1
1
1
−
sin x()
1.2
1
−
x
Рис. 1.9. Использование стандартной функции root
для решения нелинейного уравнения
sin(x) = 0
Если уравнение неразрешимо, то при попытке найти его ко-
рень будет выдано сообщение об ошибке. Кроме того, к ошибке
или выдаче неправильного корня может привести и попытка
применить метод секущих в области локального минимума или
максимума f(x). В этом случае секущая может иметь направле-
ние близкое к горизонтальному, и выводить точку следующего
приближения далеко от предполагаемого положения корня. Для
решения таких уравнений лучше применять встроенную функ-
цию Minerr. Аналогичные проблемы могут возникнуть, если на-
чальное приближение выбрано слишком далеко от настоящего
решения или f(x) имеет особенности типа бесконечности.
Иногда удобнее задавать не начальное приближение к кор-
ню, а интервал [a, b], внутри которого корень заведомо находит-
ся. В этом случае следует использовать функцию root с че-
тырьмя аргументами, а начальное значение x присваивать не
нужно. Поиск корня осуществляется в промежутке между a и b.
x root sin x() x, 1
−, 1,
(
)
:=
x 0
=
sin x() 0
=
При этом явный вид функ-
ции f(x) может быть опреде-
лен непосредственно в теле
функции root. На рис. 1.10
приведен листинг програм-
мы с использованием этого
варианта функции root.
Рис. 1.10. Поиск корня алгебраического
уравнения в заданном интервале
18
Когда функция root имеет четыре аргумента, следует пом-
нить о двух ее особенностях:
–
внутри интервала [a, b] не должно находиться более одного
корня, иначе будет найден один из них, заранее неизвестно
какой именно;
–
значения f(a) и f(b) должны иметь разный знак, иначе будет
выдано сообщение об ошибке.
Если уравнение не имеет действительных корней, но имеет
мнимые, то их также можно найти. На рис. 1.11 приведен при-
мер, в котором уравнение x
2
+ 1 = 0, имеющее два чисто мнимых
корня, решается два раза с разными начальными значениями.
x 0.5
:=
root x
2
1
+
x,
(
)
i−
=
x
0.5
−
:=
root x
2
1
+
x,
(
)
i
=
Для решения этого уравнения
второй вид функции root (с
четырьмя аргументами) не-
применим, поскольку f(x) яв-
ляется положительно опреде-
ленной и указать интервал, на
границах которого она имела
бы разный знак, невозможно.
Рис. 1.11. Поиск мнимого корня
Отметим, что f(x) может быть функцией не одного, а любо-
го количества аргументов. Эта возможность проиллюстрирована
f x y,
(
)
x
2
y
2
−
3
+
:=
x 1
:=
y 0
:=
root f x y,
(
) x,
(
)
1.732i
−
=
root f x y,
(
) y,
(
) 2
=
на рис. 1.12 на примере функ-
ции двух переменных f(x, y) =
= x
2
– y
2
+ 3. В самой функции
root необходимо определить,
относительно какого из аргу-
ментов следует решить уравне-
ние. Затем уравнение f(x, 0) = 0
решается относительно пере-
менной x, а потом
другое
уравнение – f(1, y) = 0 относи-
тельно переменной y.
Рис. 1.12. Поиск корня уравнения,
заданного функцией
двух переменных
19
При численном решении уравнений относительно одной
из переменных необходимо предварительно определить значе-
ния остальных переменных. Иначе попытка вычисления урав-
нения приведет к появлению ошибки «This variable or
function is not defined above», в данном случае говорящей о
том, что другая переменная ранее не определена. Конечно,
можно указать значения других переменных непосредственно
внутри функции root.
Если функция f(x)
–
полином, то все его корни можно оп-
ределить, используя встроенную функцию polyroots(v), где
v – вектор, составленный из коэффициентов полинома. По-
скольку полином N-й степени имеет ровно N корней (некото-
рые из них могут быть кратными), вектор v должен состоять
из
N+1
элемента. Результатом действия функции polyroots
является вектор, составленный из N корней рассматриваемого
полинома. На рис. 1.13 приведен пример решения уравнения
f(x) = (x – 13) (x – 1)
3
= x
4
– 6x
3
+12x
2
– 10x+3 = 0.
v
3
10
−
12
6
− 1
(
)
T
:=
polyroots v()
0.992
1.004 7.177i 10
3
−
×
+
1.004 7.177i 10
3
−
×
−
3
⎛
⎜
⎜
⎜
⎜
⎜
⎝
⎞
⎟
⎟
⎟
⎟
⎟
⎠
=
Коэффициенты
поли-
нома записаны в виде
вектора в первой строке
примера. Первым в век-
торе должен идти сво-
бодный член полинома,
вторым – коэффициент
при x
1
и т.д. Последним,
N + 1, элементом вектора
Рис. 1.13. Поиск корня полинома
должен быть коэффициент при старшей степени
x
N
. Во второй
строке показано действие функции polyroots. При этом числен-
ный метод вместо двух действительных единичных корней вы-
дает одинаковые мнимые числа. Однако малая мнимая часть
этих корней находится в пределах погрешности, определяемой
константой TOL, и не должна вводить пользователей в заблуж-
дение. Необходимо помнить, что корни полинома могут быть
комплексными и ошибка вычислений может сказываться как на
действительной, так и на комплексной части искомого корня.
20
В следующем примере, представленном на рис. 1.14, по-
казано вычисление трех действительных корней полинома
f(x) = 6 – 7x + x
3
с понижением порядка полинома. Здесь исполь-
зуется вариант функции root с двумя аргументами. Приведен-
ный на рисунке график функции f(x) показывает, что уравнение
имеет три действительных корня. Задавая начальное приближе-
ние z = –2, находим один из корней полинома: x
1
= –3. Затем ис-
ходный полином делится на (z – x
1
) и отыскивается второй ко-
рень x
2
= 1. Далее функция root еще раз вызывается для нахож-
дения корня полинома первого порядка, получаемого делением
исходного полинома на (z – x
1
) и (z – x
2
). Для каждого из най-
денных корней производится проверка – вычисляется невязка
уравнения.
f x
3
( )
1.904 10
4
−
×
=
x
3
2
=
x
3
root
f z()
z x
1
−
( )
z x
2
−
( )
⋅
z,
⎡
⎢
⎣
⎤
⎥
⎦
:=
z 3
:=
f x
2
( )
6.078 10
4
−
×
=
x
2
1
=
x
2
root
f z()
z x
1
−
z,
⎛
⎜
⎝
⎞
⎟
⎠
:=
z 0.5
:=
f x
1
( )
2.463 10
5
−
×
=
x
1
3−
=
x
1
root f z() z,
(
)
:=
z
2−
:=
f y()
a
1
a
2
y⋅
+
a
3
y
2
⋅
+
a
4
y
3
⋅
+
:=
a
4
1
:=
a
3
0
:=
a
2
7−
:=
a
1
6
:=
2
0
2
10
0
10
f y()
Рис. 1.14. Поиск корня полинома с понижением порядка
21
2. МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
2.1. Общие вопросы
По оценкам современной научной литературы, около 75 %
всех вычислительных задач приводят к решению систем линей-
ных алгебраических уравнений (СЛАУ). Иногда СЛАУ получа-
ется непосредственно как математическая модель какого-либо
процесса, иногда является приближением (аппроксимацией)
дифференциальных/интегральных уравнений, описывающих
физический процесс. Поэтому для решения задачи надо уметь
быстро и качественно решать СЛАУ. Конечно, существует мно-
го методов и современных пакетов прикладных программ для
решения СЛАУ, но для того, чтобы их успешно применять, не-
обходимо разбираться в основах построения методов и алгорит-
мов, иметь представление о недостатках и преимуществах ис-
пользуемых методов.
Все методы решения линейных алгебраических задач (на-
ряду с задачей на решение СЛАУ – это и вычисление определи-
телей, и обращение матриц, и задачи на собственные значения)
можно разбить на два класса: прямые (точные) и итерационные
(приближенные).
Прямые методы позволяют получить решение за конечное
число арифметических операций. Если операции реализуются
точно, то и решение будет точным (поэтому прямые методы на-
зываются еще точными методами). В итерационных методах
решением является предел некоторой бесконечной последова-
тельности единообразных действий.
Будем решать лишь такие СЛАУ, у которых число уравне-
ний совпадает с числом неизвестных, причем будем предпола-
гать наличие единственного решения.
Задана СЛАУ размерности m:
⎪
⎪
⎩
⎪
⎪
⎨
⎧
=
⋅
+
+
⋅
+
⋅
=
⋅
+
+
⋅
+
⋅
=
⋅
+
+
⋅
+
⋅
m
m
mm
m
m
m
m
m
m
f
x
a
x
a
x
a
f
x
a
x
a
x
a
f
x
a
x
a
x
a
...
...
...
...
2
2
1
1
2
2
2
22
1
21
1
1
2
12
1
11
,
(2.1)
22
которую можно записать также в матричном виде:
f
x
r
r
=
A
,
(2.1′)
где
⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎛
=
mm
m
m
m
m
a
a
a
a
a
a
a
a
a
...
...
...
A
2
1
2
22
21
1
12
11
L
L
L
L
,
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎝
⎛
=
−
m
m
f
f
...
f
f
f
1
2
1
r
,
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎝
⎛
=
−
m
m
x
x
...
x
x
x
1
2
1
r
,
A – матрица системы,
f
r
–
вектор правых частей, x
r
–
вектор не-
известных. Система имеет решение, если det A ≠ 0. По опреде-
лению решения, подставив вектор
x
r
в СЛАУ, получим m тож-
дественных уравнений.
Эффективность способов решения системы (2.1) во многом
зависит от структуры и свойств матрицы A: размерности, обу-
словленности, симметричности, заполненности (т.е. соотноше-
ния между числом ненулевых и нулевых элементов) и др.
2.2. Точные методы решения систем линейных
алгебраических уравнений
Метод Крамера
При небольшой размерности системы m (m = 2, 3) на прак-
тике часто используют формулы Крамера решения СЛАУ:
A
A
det
det
i
i
x =
(i = 1, 2, …, m).
Эти формулы позволяют находить неизвестные в виде дро-
бей, знаменателем которых является определитель матрицы сис-
темы, а числителем – определители матрицы A
i
, полученные из
A заменой столбца коэффициентов при вычисляемом неизвест-
ном столбцом свободных членов.
Размерность системы (т.е. число m) является главным фак-
тором, из-за которого формулы Крамера не могут быть исполь-
зованы для численного решения СЛАУ большого порядка. При
непосредственном раскрытии определителей решение системы с
m неизвестными требует порядка m!m арифметических опера-
23
ций. Таким образом, для решения системы, например, из
m = 100 уравнений потребуется совершить 10
158
операций, что
не под силу даже самым мощным современным ЭВМ.
Метод обратной матрицы
Если det A ≠ 0, то существует обратная матрица A
–1
. По оп-
ределению обратной матрицы, это такая матрица, что
A A
–1
= A
–1
A = E, где E – единичная матрица:
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎝
⎛
=
1
...
0
0
0
...
1
0
0
...
0
1
K
K
K
K
E
.
Если обратная матрица известна, то, умножая на нее СЛАУ
слева, получим:
.
,
,
1
1
1
1
f
x
f
x
f
x
r
r
r
r
r
r
−
−
−
−
=
=
=
A
A
E
A
A
A
Следо-
вательно, решение СЛАУ свелось к умножению известной об-
ратной матрицы на вектор правых частей. Таким образом, зада-
ча решения СЛАУ и задача нахождения обратной матрицы свя-
заны между собой, поэтому часто решение СЛАУ называют за-
дачей обращения матрицы.
Метод обратной матрицы можно использовать для решения
систем небольших размерностей. На рис. 2.1 показано решение
системы с заданной матрицей A и вектором правых частей b с
применением стандартных функций нахождения обратной мат-
рицы и умножения матрицы на вектор.
A
8
3
3
4
5
2−
2
1
10
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
:=
b
10
5
4
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
:=
x A
1
−
b⋅
:=
x =
Рис. 2.1. Решение СЛАУ методом обратной матрицы
Метод Гаусса
Наиболее известным и популярным точным способом ре-
шения линейных систем вида (2.1) является метод Гаусса. Этот
метод заключается в последовательном исключении неизвест-
ных. Пусть в системе уравнений
24
⎪
⎪
⎩
⎪
⎪
⎨
⎧
=
⋅
+
+
⋅
+
⋅
=
⋅
+
+
⋅
+
⋅
=
⋅
+
+
⋅
+
⋅
)
0
(
)
0
(
2
)
0
(
2
1
)
0
(
1
)
0
(
2
)
0
(
2
2
)
0
(
22
1
)
0
(
21
)
0
(
1
)
0
(
1
2
)
0
(
12
1
)
0
(
11
...
...
...
...
m
m
mm
m
m
m
m
m
m
f
x
a
x
a
x
a
f
x
a
x
a
x
a
f
x
a
x
a
x
a
первый элемент
0
)0
(
11
≠
a
. Назовем его ведущим элементом пер-
вой строки. Поделим все элементы этой строки на
)0
(
11
a
и исклю-
чим x
1
из всех последующих строк, начиная со второй, путем
вычитания первой (преобразованной), умноженной на коэффи-
циент при
1
x в соответствующей строке.
Получим
⎪
⎪
⎩
⎪
⎪
⎨
⎧
=
⋅
+
+
⋅
+
⋅
=
⋅
+
+
⋅
+
⋅
=
⋅
+
+
⋅
+
⋅
+
)1
(
)1
(
3
)1
(
3
2
)1
(
2
)1
(
2
)1
(
2
3
)1
(
23
2
)1
(
22
)1
(
1
)1
(
1
3
)1
(
13
2
)1
(
12
1
...
...
...
...
m
m
mm
m
m
m
m
m
m
f
x
a
x
a
x
a
f
x
a
x
a
x
a
f
x
a
x
a
x
a
x
.
Если
0
)1
(
22
≠
a
, то, продолжая аналогичное исключение, при-
ходим к системе уравнений с верхней треугольной матрицей
⎪
⎪
⎪
⎩
⎪
⎪
⎪
⎨
⎧
=
=
⋅
+
+
=
⋅
+
+
⋅
+
=
⋅
+
+
⋅
+
⋅
+
)
(
)3
(
3
)3
(
3
3
)2
(
2
)2
(
2
3
)2
(
23
2
)1
(
1
)1
(
1
3
)1
(
13
2
)1
(
12
1
...
...
...
m
m
m
m
m
m
m
m
m
f
x
f
x
a
x
f
x
a
x
a
x
f
x
a
x
a
x
a
x
K
K
K
.
Из нее в обратном порядке находим все значения x
i
:
⎪
⎪
⎩
⎪
⎪
⎨
⎧
⋅
−
−
⋅
−
⋅
−
=
⋅
−
=
=
−
−
−
−
−
m
m
m
m
m
m
m
m
m
m
m
m
x
a
x
a
x
a
f
x
x
a
f
x
f
x
)1
(
1
3
)1
(
13
2
)1
(
12
)1
(
1
1
)1
(
1
)1
(
1
1
)
(
...
K
K
K
.
25
Процесс приведения к системе с треугольной матрицей на-
зывается прямым ходом, а нахождения неизвестных – обратным.
Если один из ведущих элементов равен нулю, изложенный ал-
горитм метода Гаусса неприменим. Кроме того, если какие-либо
ведущие элементы малы, то это приводит к увеличению ошибок
округления и ухудшению точности счета. Поэтому обычно ис-
пользуется другой вариант метода Гаусса – схема Гаусса с вы-
бором главного элемента. Путем перестановки строк, а также
столбцов с соответствующей перенумерацией коэффициентов и
неизвестных добиваются выполнения условия:
)0
(
)0
(
ij
ii
a
a
≥
, i, j = 1, 2, …, m,
т.е. осуществляется выбор первого главного элемента. Разделив
первую строку на главный элемент, как и прежде, исключают x
1
из остальных уравнений. Затем для оставшихся столбцов и
строк выбирают второй главный элемент и т.д.
Метод прогонки
Часто возникает необходимость в решении СЛАУ, матрицы
которых являются слабо заполненными, т.е. содержат много ну-
левых элементов. В то же время эти матрицы имеют определен-
ную структуру. Среди таких систем выделим системы с матри-
цами ленточной структуры, в которых ненулевые элементы рас-
полагаются на главной диагонали и на нескольких побочных
диагоналях. Для решения систем с ленточными матрицами ко-
эффициентов вместо метода Гаусса можно использовать более
эффективные методы.
Рассмотрим наиболее простой случай: систему с трехдиаго-
нальной матрицей коэффициентов, к которой сводится решение
ряда численных задач (сплайн-интерполяция таблично заданной
функции, дискретизация краевых задач для дифференциальных
уравнений методами конечных разностей и др.). Тогда СЛАУ
можно записать в упрощенном виде:
i
i
i
i
i
i
i
f
x
b
x
c
x
a
=
+
−
+
−
1
1
,
(2.2)
где i = 1, 2, …, m. Схема (2.2) имеет трехдиагональную структу-
ру, что хорошо видно из следующего, эквивалентного (2.2), век-
торно-матричного представления:
26
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
⋅
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
−
−
−
−
−
−
−
−
−
m
m
m
m
m
m
m
m
m
f
f
f
f
f
x
x
x
x
x
c
a
b
c
a
b
c
a
b
c
a
b
c
1
3
2
1
1
3
2
1
1
1
1
3
3
3
2
2
2
1
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
M
M
K
K
K
K
K
K
K
K
K
K
K
K
K
.
Если при этом выполняется условие
i
i
i
a
b
c
+
≥
, то говорят, что
матрица этой системы имеет диагональное преобладание.
Предположим, что существуют такие наборы чисел
α
i
и
β
i
(i = 1, 2, …, m–1), при которых
1
1
1
+
+
+
+
=
i
i
i
i
x
x
β
α
.
(2.3)
Уменьшим в (2.3) индекс на единицу:
i
i
i
i
x
x
β
α
+
=
−1
и подставим полученное выражение в (2.2)
i
i
i
i
i
i
i
i
i
i
f
x
b
x
c
a
x
a
=
+
−
+
+1
β
α
,
откуда
i
i
i
i
i
i
i
i
i
i
i
i
a
c
f
a
x
a
c
b
x
α
β
α
−
−
+
−
=
+1
. Данное равенство в точности
совпадает с (2.3), если при всех i = 1, 2, …, m–1 выполняются
рекуррентные соотношения
i
i
i
i
i
a
c
b
α
α
−
=
+1
,
i
i
i
i
i
i
i
a
c
f
a
α
β
β
−
−
=
+1
.
(2.4)
Процесс вычисления
α
i
и
β
i
можно начать со значений
1
1
1
c
b
=
α
,
1
1
1
c
f
−
=
β
.
(2.4′)
Далее продолжаем по формулам (2.4) последовательно при
i = 1, 2, …, m–1. При i = m из (2.2) имеем
m
m
m
m
m
f
x
c
x
a
=
−
−1
.
(2.5)
В то же время при i = m–1 из (2.3) получаем
m
m
m
m
x
x
β
α
+
=
−1
.
(2.6)
Подставляя выражение для x
m–1
из (2.6) в (2.5) и решая получен-
ное выражение относительно x
m
, получаем:
27
m
m
m
m
m
m
m
a
c
f
a
x
α
β
−
−
=
,
(2.7)
где
α
m
и
β
m
известны из предыдущего шага. Далее по формулам
(2.3) последовательно находятся x
m–1
, x
m–2
, …, x
1
.
Данный способ решения системы уравнений вида (2.2) назы-
вается методом прогонки. Исходя из вышеизложенного, можно
выделить три этапа метода прогонки. На первом этапе по форму-
лам (2.4′), (2.4) определяются так называемые прогоночные ко-
эффициенты
α
i
и
β
i
при i = 2, 3, …, m (прямая прогонка). Затем
по формуле (2.7) находится x
m
. На последнем этапе по формуле
(2.3) определяются x
i
при i = m–1, m–2, …, 1 (обратная прогонка).
Для успешного применения метода прогонки нужно, чтобы
в процессе вычислений не возникало ситуаций с делением на
нуль, а при больших размерностях систем не должно быть бы-
строго роста погрешностей округления.
Будем называть прогонку корректной, если знаменатели
прогоночных коэффициентов (2.4) не обращаются в нуль, и ус-
тойчивой, если |
α
i
| < 1 при i = 1, 2, 3, …, m.
Теорема 2.1. Пусть коэффициенты a
i
, b
i
уравнения (2.2) при
i = 2, 3, …, m–1 отличны от нуля и пусть
i
i
i
a
b
c
+
≥
при i = 1, 2, 3, …, m.
Тогда прогонка (2.3)–(2.7) корректна и устойчива.
Условия этой теоремы, которые во многих приложениях
выполняются автоматически, являются достаточными условия-
ми корректности и устойчивости прогонки. Если эти условия не
выполняются, то можно аналогично схеме Гаусса организовать
выбор главного элемента.
2.3. Итерационные методы решения линейных
алгебраических систем
Решение СЛАУ методом простых итераций
Система (2.1′) может быть преобразована к эквивалентной
ей системе вида:
β
r
r
r
+
= x
x α
,
(2.8)
28
где x
r
– искомый вектор, а
α
и
β
r
– новые матрица и вектор. Бу-
дем решать (2.8) методом последовательных приближений. За-
дадим нулевое приближение:
( )
b
x
r
r
=
0
, тогда
( )
1
x
r
определяется
рекуррентным равенством
( )
( )
β
r
r
r
+
=
0
1
x
x
α
,
далее находим
( )
2
x
r
( )
( )
β
r
r
r
+
=
1
2
x
x
α
.
Для k-й итерации получаем
(
)
( )
...
,2
,1,
0
,
1
=
+
=
+
k
x
x
k
k
β
r
r
r
α
(2.9)
Такой итерационный процесс будем называть методом простых
итераций (МПИ). Изучим вопрос о сходимости этого процесса,
т.е. определим, какие нужно предъявить требования к виду мат-
рицы
α
, чтобы
( )
*
lim
x
x
k
k
r
r
=
∞
→
.
Теорема 2.2. Необходимым и достаточным условием схо-
димости МПИ (2.9) при любом начальном векторе
( )
0
x
r
к
решению системы (2.8) является требование, чтобы норма
матрицы
α
была меньше 1:
1
<
α
.
(2.10)
Это требование равносильно условию малости элементов
матрицы
α
по абсолютной величине, т.к. в качестве нормы мат-
рицы используют максимальное значение из сумм модулей эле-
ментов строк этой матрицы
∑
=
≤
≤
=
n
j
ij
n
i
a
1
1
max
α
.
Важной проблемой является вопрос о способе остановки
итерационного процесса при достижении точности. Наиболее
простой способ – это сравнение между собой соответствующих
неизвестных на двух соседних итерациях (k+1) и (k). Если мак-
симальная из всех разностей становится меньше заданной точ-
ности
ε
, то итерационный процесс останавливается
29
ε
<
−
+
≤
≤
1
1
max
k
i
k
i
n
i
x
x
.
Можно применить способ, связанный с вычислением вектора
невязки r
r
:
i
n
j
k
ij
i
b
x
a
r
j
−
=
∑
=1
, показывающий, насколько полу-
ченное приближение
k
x
r
отличается от точного решения. Затем
вычисляется норма вектора невязки
|
|
max
1
i
n
i
r
r
≤
≤
=
r
. Если она ма-
ла, т.е.
ε
<
r
r
, то итерационный процесс останавливается.
Рассмотрим несколько способов построения МПИ.
Метод Якоби
Предположим, что диагональные элементы матрицы A ис-
ходной системы (2.1′) не равны 0 (a
ii
≠ 0, i = 1, 2, …, n). Разре-
шим первое уравнение системы (2.1) относительно x
1
, второе
относительно x
2
и т.д. Получим систему в виде (2.8):
⎪
⎪
⎩
⎪
⎪
⎨
⎧
+
+
+
+
=
+
+
+
+
=
+
+
+
+
=
−
−
n
n
n
n
n
n
n
n
n
n
n
x
x
x
x
x
x
x
x
x
x
x
x
β
α
α
α
β
α
α
α
β
α
α
α
1
1
,
2
2
1
1
2
2
3
23
1
21
2
1
1
3
13
2
12
1
...
..
..........
..........
..........
..........
..........
...
...
,
где
n
i
a
b
j
i
j
i
a
a
ii
i
i
ii
ij
ij
...,
,2
,1
,
,
,0
,
=
=
⎪
⎩
⎪
⎨
⎧
=
≠
−
=
β
α
.
Метод, основанный на таком приведении СЛАУ к виду
(2.8), называют методом Якоби. Теперь, задав нулевое прибли-
жение, по рекуррентным соотношениям (2.9) можем выполнять
итерационный процесс. Сформулированное выше условие схо-
димости в методе Якоби равносильно условию диагонального
преобладания:
.
...,
,2
,1
,
1
n
i
a
a
n
j
i
j
ij
ii
=
>
∑
≠
=
30
α
0
3
−
5
3
−
10
1
−
2
0
1
5
1
−
4
1
−
5
0
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠
:=
β
10
8
1
4
10
⎛
⎜
⎜
⎜
⎜
⎜
⎝
⎞
⎟
⎟
⎟
⎟
⎟
⎠
:=
yakobi ε
( )
err
100
←
x
β
←
x1
α x
⋅
β
+
←
err
x1 x
−
←
x
x1
←
err
ε
>
while
x
:=
yakobi 0.001
(
) =
Действительно, если это
условие выполняется, то и
суммы модулей элементов
строк матрицы α меньше 1.
На рис. 2.2 приведен
пример решения в MathCAD
методом Якоби системы с
матрицей A и вектором пра-
вых частей b:
.
4
5
10
,
10
2
3
1
5
3
2
4
8
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎝
⎛
=
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎝
⎛
−
=
b
A
Легко убедиться, что для
исходной матрицы A вы-
полняются условия диаго-
нального преобладания, а
для матрицы α – условия
(2.11), что
обеспечивает
сходимость итерационного
процесса.
Рис. 2.2. Решение СЛАУ
методом Якоби
Метод Зейделя
Под методом Зейделя обычно понимается такое видоизме-
нение МПИ (2.10) решения СЛАУ (2.8), в котором для подсчета
i-й компоненты (k+1)-го приближения к искомому вектору
*
x
r
используются уже вычисленные на этом, т.е. (k+1)-м шаге, но-
вые значения первых i–1 компонент. Это означает, что если сис-
тема (2.8) тем или иным способом сведена (например, с помо-
щью метода Якоби) к системе (2.9) с матрицей коэффициентов
α и вектором свободных членов
,β
r
то ее приближение к реше-
нию по методу Зейделя определяется системой равенств
(
)
( )
( )
( )
( )
(
)
(
)
( )
( )
( )
(
)
(
)
(
)
(
)
( )
⎪
⎪
⎩
⎪
⎪
⎨
⎧
+
+
+
+
+
=
+
+
+
+
+
=
+
+
+
+
+
=
+
−
−
+
+
+
+
+
+
3
,
1
1
1
,
1
2
2
1
1
1
1
2
2
3
23
2
22
1
1
21
1
2
1
1
3
13
2
12
1
11
1
1
...
..........
..........
..........
..........
..........
..........
..........
..........
...
...
β
α
α
α
α
β
α
α
α
α
β
α
α
α
α
k
n
n
n
k
n
n
n
k
n
k
n
k
n
k
n
n
k
k
k
k
k
n
n
k
k
k
k
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
. (2.11)
31
С точки зрения компьютерной реализации МПИ, использо-
вание метода Зейделя означает, что элементы массива x
r
будут
постепенно замещаться новыми элементами. В связи с такой ин-
терпретацией метод Зейделя иногда называют методом после-
довательных смещений.
Метод релаксации
Иногда исходную систему (2.8) не удается привести к виду
(2.9), выполняя при этом условие сходимости (2.11). В этом слу-
чае можно воспользоваться методом релаксации, который осно-
вывается на соотношении
( )
( )
( )
b
x
x
x
k
k
k
r
r
r
r
+
−
=
−
+
A
τ
1
, откуда
( )
( )
( )
(
)
b
x
x
x
k
k
k
r
r
r
r
−
−
=
+
A
τ
1
, где
τ
– параметр релаксации. Скаляр-
ные формулы метода релаксации имеют следующий вид:
(
)
( )
( )
( )
( )
(
)
(
)
( )
( )
( )
( )
(
)
(
)
( )
( )
( )
( )
(
)
⎪
⎪
⎩
⎪
⎪
⎨
⎧
−
+
+
+
−
=
−
+
+
+
−
=
−
+
+
+
−
=
+
+
+
n
k
n
nn
k
n
k
n
k
k
n
k
n
n
k
k
k
k
k
n
n
k
k
k
k
b
x
a
x
a
x
a
x
x
b
x
a
x
a
x
a
x
x
b
x
a
x
a
x
a
x
x
...
........
..........
..........
..........
..........
..........
..........
..........
...
...
2
2
1
1
1
1
2
2
2
22
1
21
2
1
2
1
1
2
12
1
11
1
1
1
τ
τ
τ
.
(2.12)
Раскрыв скобки, можно привести (2.12) к виду (2.10), где ко-
эффициенты матрицы α и вектор свободных членов
β
r
будут
иметь вид:
n
...,
,
,
i
,
b
,
j
i
,
a
j
i,
a
i
i
ij
ij
ij
2
1
1
=
=
⎩
⎨
⎧
≠
−
=
−
=
τ
β
τ
τ
α
. Подбором
параметра
τ
можно добиться сходимости метода релаксации.
2.4. Стандартные функции пакета MathCAD
В MathCAD СЛАУ можно решить как в развернутой форме
(2.1), так и в более компактной форме (2.2). Для первого способа
следует использовать вычислительный блок Given/Find, со-
стоящий из трех последовательных частей:
–
Given – ключевое слово;
–
система, записанная логическими операторами в виде ра-
венств и, возможно, неравенств;
–
Find(x
1
, … , x
M
) – встроенная функция для решения систе-
мы относительно переменных x
1
, … , x
M
.
32
Перед вызовом вычислительного блока всем неизвестным
присвоены начальные значения. Они могут быть произвольными,
т.к. решение СЛАУ с невырожденной матрицей единственно.
Для второй формы записи системы используют встроенную
функцию lsolve:
lsolve(A, b) – решение системы линейных уравнений;
A – матрица коэффициентов системы;
b – вектор правых частей.
На рис. 2.3–2.5 приведены примеры решения СЛАУ с по-
мощью стандартных функций MathCAD.
A
1
0.7
3
5
12
0
2
5
4
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
:=
b
1
2.9
3.1
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
:=
lsolve A b,
(
)
0.186
−
0.129
−
0.915
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
=
Рис. 2.3. Решение СЛАУ
с помощью вычислительного
блока Given/Find
Рис. 2.4. Решение СЛАУ
в форме (2.2)
A
1
0.7
3
5
12
0
2
5
4
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
:=
b
1
2.9
3.1
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
:=
lsolve A b,
(
)
0.186
−
0.129
−
0.915
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
=
Рис. 2.5. Символьное решение СЛАУ
x 0
:=
y 0
:=
z
0
:=
Given
1 x 5y
2 z 1
0.7 x 12y 5 z 2.9
3 x 0y 4 z 3.1
Find x y z
(
)
0.186
−
0.129
−
0.915
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎟
⎠
=
33
3. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ
Слово «интерполяция» в переводе с латыни означает «меж-
ду точками». Задачи интерполяции часто возникают в инженер-
ных и других практических приложениях. Допустим, что в ре-
зультате экспериментальных измерений получена таблица зна-
чений некоторой функции. Требуется найти промежуточные
значения этой функции, а также производные, определяющие
скорость ее изменения. Это так называемая задача о восстанов-
лении функции. Кроме того, при проведении расчетов сложные
функции удобно заменять алгебраическими многочленами или
другими элементарными функциями, которые достаточно про-
сто вычисляются (задача о приближении функции). Интерполя-
цию используют для приближенного вычисления интегралов
(построение квадратурных формул). Из математического ана-
лиза известны, например, многочлены (ряды) Тейлора, которые
применяют для вычисления значений гладких (т.е. достаточное
число раз дифференцируемых) функций. В точных науках часто
используют разложение функций в тригонометрические ряды.
Каждый метод имеет свою погрешность, определяемую тем, на-
сколько различаются значения исходной и интерполирующей
функций. Существуют ли другие способы интерполяции и при-
ближения функций? Когда и какой способ лучше применять?
Какова точность (погрешность) используемых методов интер-
поляции? Об этом мы узнаем, изучив следующую тему и вы-
полнив соответствующую лабораторную работу.
3.1. Постановка задачи интерполяции
На интервале [a, b] задана система узлов интерполяции
x
i
, i = 0, 1,..., N, a ≤ x
i
≤ b, и значения неизвестной функции в
этих узлах f
i
, i = 0, 1,..., N. Могут быть поставлены задачи:
1. Найти функцию F(x), принимающую в точках
x
i
заданные
значения: F(x
i
) = f
i
, i = 0, 1,…, N (условия интерполяции).
2. Для заданного значения z ∈ [a, b] найти F(z).
3. Для заданного значения z ∈ [a, b] найти F′(z).
Задача имеет много решений: через заданные точки (x
i
, f
i
),
i = 0, 1,..., N,
можно провести бесконечно много кривых, каждая
из которых будет графиком функции, для которой выполнены
34
все условия интерполяции. Если известна исходная функция
g(x), то можно оценить погрешность метода в произвольной
точке z ∈ [a, b]: r(z) = |g(z) – F(z)|. Кроме того, можно оценивать
равномерную r
1
и среднеквадратичную r
2
погрешности:
)
(
max
]
,
[
1
z
r
r
b
a
z∈
=
,
(
)
dz
)z
(r
r
b
a
∫
=
2
2
.
Нас будет интересовать поведение погрешности метода при
увеличении числа узлов интерполяции. Будем говорить, что ме-
тод сходится, если при N→
∞
погрешность r → 0.
Все методы интерполяции можно разделить на локальные и
глобальные. В случае локальной интерполяции на каждом ин-
тервале [x
i–1
, x
i
] строится своя (локальная) функция. В случае
глобальной интерполяции отыскивается одна (глобальная)
функция на всем интервале [a, b]. Далее приведены примеры
различных способов интерполяции.
3.2. Локальная интерполяция
Кусочно-постоянная интерполяция. На каждом локальном
отрезке [x
i–1
, x
i
], i = 1, 2,…, N, интерполирующая функция явля-
ется постоянной и равна левому: F
i
(z) = f
i
или правому: F
i
(z) = f
i
значению. Легко понять, что условия интерполяции выполняют-
ся. Построенная функция разрывная (рис. 3.1а), что ограничива-
ет ее применение. Кроме того, в случае малого числа точек та-
кая интерполяция дает большую погрешность.
Кусочно-линейная интерполяция. На каждом интервале
[x
i–1
, x
i
] функция является линейной F
i
(z) = k
i
z+l
i
. Значения ко-
эффициентов находятся из выполнения условий интерполяции на
концах отрезка: F
i
(x
i–1
) = f
i–1
, F
i
(x
i
) = f
i
. Получаем систему урав-
нений: k
i
x
i–1
+l
i
= f
i–1
, k
i
x
i
+l
i
= f
i
,, откуда находим
,
1
1
−
−
−
−
=
i
i
i
i
i
x
x
f
f
k
i
i
i
i
x
k
f
l
−
=
. Итоговая функция будет непрерывной, но произ-
водная будет разрывной в каждом узле интерполяции. Погреш-
ность такой интерполяции будет меньше, чем в предыдущем
случае. Иллюстрация кусочно-линейной интерполяции приве-
дена на рис. 3.1б.
35
x
y
x
N
x
i
f
0
f
N
f
i
x
0
а
x
y
x
N
x
i–1
f
0
f
N
x
0
б
x
i
f
i–1
f
i
Рис. 3.1. Левая кусочно-постоянная (а)
и кусочно-линейная (б) интерполяции
Кусочно-параболическая интерполяция. На каждом i-м ин-
тервале [x
i–1
, x
i
], i = 1, 2,…, N, функция является параболой вида
F
i
(z) = a
i
+ b
i
(z – x
i
) + c
i
/2(z – x
i
)
2
. Значения коэффициентов a
i
, b
i
,
c
i
находятся из условий интерполяции, непрерывности функции
и ее первой производной в узлах интерполяции.
Кубический интерполяционный сплайн. Название метода про-
исходит от английского слова «spline». Так называлась гибкая
линейка, использовавшаяся корабельными инженерами вместо
лекал. Форма этого универсального лекала на каждом отрезке
описывается кубической параболой. Сплайны широко применя-
ются в инженерных приложениях, в частности в компьютерной
графике. Итак, на каждом i-м интервале [x
i–1
, x
i
], i = 1, 2,…, N,
решение будем искать в виде полинома третьей степени:
S
i
(x) = a
i
+ b
i
(x – x
i
) + c
i
(x – x
i
)
2
/2 + d
i
(x – x
i
)
3
/6.
Неизвестные коэффициенты a
i
, b
i
, c
i
, d
i
, i = 1, 2,..., N, находим:
– из условий интерполяции: S
i
(x
i
) = f
i
, i = 1, 2,..., N, S
1
(x
0
) = f
0
;
– непрерывности функции S
i
(x
i–1
) = S
i–1
(x
i–1
), i = 2, 3,..., N;
– непрерывности первой и второй производной:
S
′
i
(x
i–1
) = S
′
i–1
(x
i–1
), S
′′
i
(x
i–1
) = S
′′
i–1
(x
i–1
), i = 2, 3,..., N.
Для определения 4N неизвестных получаем систему 4N–2 урав-
нений:
a
i
= f
i
, i = 1, 2,..., N,
b
i
h
i
– c
i
h
i
2
/2 + d
i
h
i
3
/6 = f
i
– f
i–1
, i = 1, 2,..., N,
b
i
– b
i–1
= c
i
h
i
– d
i
h
i
2
/2, i = 2, 3,..., N,
36
d
i
h
i
= c
i
– c
i–1
, i = 2, 3,..., N,
где h
i
= x
i
– x
i–1
. Недостающие два уравнения выводятся из до-
полнительных условий, например S"(a) = S"(b) = 0. Из системы
можно исключить неизвестные b
i
, d
i
, получив систему N+1 ли-
нейных уравнений (СЛАУ) для определения коэффициентов c
i
:
c
0
= 0, c
N
= 0,
h
i
c
i–1
+2(h
i
+h
i+1
)c
i
+h
i+1
c
i+1
= 6
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
−
−
−
−
+
1
1
1
i
i
i
i
i
i
h
f
f
h
f
f
, i = 1, 2,…, N–1.
Данная CЛАУ имеет трехдиагональную матрицу и решается ме-
тодом прогонки. После этого вычисляются коэффициенты b
i
, d
i
:
i
i
i
i
i
i
i
i
i
i
i
i
h
f
f
h
d
h
c
,b
h
c
c
d
1
2
1
6
2
−
−
−
+
−
=
−
=
, i = 1, 2,..., N.
Для вычисления значения S(x) в произвольной точке отрез-
ка z∈ [a, b] необходимо решить систему уравнений на коэффи-
циенты c
i
, i = 1, 2,…, N–1, затем найти все коэффициенты b
i
, d
i
.
Далее необходимо определить, на какой интервал [x
i0,
x
i0–1
] по-
падает эта точка, и, зная номер i
0
, вычислить значение сплайна и
его производных в точке z:
S(z) = a
i0
+ b
i0
(z – x
i0
) + c
i0
(z – x
i0
)
2
/2 + d
i0
(z – x
i0
)
3
/6,
S
′
(z) = b
i0
+ c
i0
(z – x
i0
) + d
i0
(z – x
i0
)
2
/2, S
′′
(z) = c
i0
+ d
i0
(z – x
i0
).
Кубический интерполяционный сплайн имеет достаточно
хорошую точность и в то же время простую и экономичную
реализацию (метод прогонки).
Для построения кубического интерполяционного сплайна в
MathCAD определена функция interp(s,x,y,t), аппроксимирующая
данные векторов x и y кубическими сплайнами. Здесь s – вектор
вторых производных от табличной функции (x, y), созданный
одной из функций cspline(x,y), pspline(x,y) или lspline(x,y);
x – вектор действительных данных аргумента, элементы которо-
го расположены в порядке возрастания; y – вектор действи-
тельных данных значений того же размера; t – значение аргу-
мента, при котором вычисляется интерполирующая функция.
Встроенные функции lspline(x,y), pspline(x,y) и cspline(x,y)
различаются способом вычисления второй производной в гра-
ничных точках. На рис. 3.2 приведен пример сплайн-аппрокси-
37
мации функции y = cos(x) на отрезке [0, 6] с шагом h = 1. Для
вычисления вторых производных использована функция
cspline. Кроме значений самого сплайна, вычислены также его
первая и вторая производные. Сравнение точных значений ис-
ходной функции, ее первой (z) и второй (u) производных в точ-
ках массива xx, расположенных между узлами интерполяции, с
вычисленными по сплайн-интерполяции значениями показыва-
ет, что, несмотря на очень грубую сетку (h = 1), значения функ-
ции и ее производных восстанавливаются достаточно хорошо.
f x()
cos x()
:=
x
0 1 2 3 4 5 6
(
)
T
:=
y
f x()
:=
s
cspline x y,
(
)
:=
A t()
interp s x, y, t,
(
)
:=
B t()
t
interp s x, y, t,
(
)
d
d
:=
C t()
2
t
interp s x, y, t,
(
)
d
d
2
:=
xx
0.5
1.5
2.5
3.5
4.5
5.5
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠
:=
A xx
( ) =
f xx
( ) =
pr_sp
h
k
B xx
k
( )
←
k 0 5
..
∈
for
h
:=
pr2_sp
g
k
C xx
k
( )
←
k 0 5
..
∈
for
g
:=
z
sin xx
( )
−
:=
u
cos xx
( )
−
:=
pr_sp =
z =
pr2_sp =
u =
Рис. 3.2. Построение кубического интерполяционного сплайна
38
3.3. Глобальная интерполяция
Будем искать интерполирующую функцию в виде полинома
(многочлена) m-й степени P
m
(x) = a
0
+ a
1
x + a
2
x
2
+ a
3
x
3
+ … + a
m
x
m
.
Какова должна быть степень многочлена, чтобы удовлетворить
всем условиям интерполяции? Допустим, что заданы две точки:
(x
0
, f
0
) и (x
1
, f
1
). Через эти точки можно провести единственную
прямую, т.е. интерполирующей функцией будет полином пер-
вой степени P
1
(x) = a
0
+ a
1
x. Через три точки можно провести
параболу P
2
(x) = a
0
+ a
1
x + a
2
x
2
и т.д. Рассуждая таким образом,
можно предположить, что искомый полином должен иметь сте-
пень N.
Чтобы доказать это, выпишем систему уравнений на коэф-
фициенты. Уравнения системы представляют собой условия ин-
терполяции при каждом x = x
i
:
( )
( )
( )
⎪
⎪
⎩
⎪
⎪
⎨
⎧
=
+
…
+
+
+
+
=
=
+
…
+
+
+
+
=
=
+
…
+
+
+
+
=
N
N
N
N
N
N
N
N
N
N
N
N
N
N
N
f
x
a
x
a
x
a
x
a
a
x
P
f
x
a
x
a
x
a
x
a
a
x
P
f
x
a
x
a
x
a
x
a
a
x
P
3
3
2
2
1
0
2
2
3
2
3
2
2
2
2
1
0
2
1
1
3
1
3
2
1
2
1
1
0
1
...
.
Данная система является линейной относительно искомых
коэффициентов a
0
, a
1
, a
2
,…, a
N
. Известно, что СЛАУ имеет ре-
шение, если ее определитель отличен от нуля. Определитель
данной системы
(
)
∏
≤
<
≤
−
=
=
Δ
N
m
k
m
k
N
N
N
N
N
x
x
x
x
...
...
x
...
x
x
...
x
0
2
2
1
1
1
1
1
носит имя Вандермонда. Из курса математического анализа из-
вестно, что он отличен от нуля, если x
k
≠ x
m
(т.е. все узлы интер-
поляции различные). Таким образом, доказано, что система
имеет решение.
Мы показали, что для нахождения коэффициентов
a
0
, a
1
, a
2
,…, a
N
надо решить СЛАУ, что является сложной зада-
39
чей. Но есть другой способ построения полинома N-й степени,
который не требует решения такой системы.
3.4. Полином Лагранжа
Решение ищем в виде
( )
( )
z
lf
z
L
i
N
i
i
n
∑
=
=
0
, где l
i
(z) – базисные
полиномы N-й степени, для которых выполняется условие:
⎩
⎨
⎧
≠
=
=
k
i
k
i
x
l
k
i
,0
,1
)
(
. Убедимся в том, что если такие полиномы по-
строены, то L
N
(x) будет удовлетворять условиям интерполяции:
( )
( )
( )
( )
( )
( )
i
i
N
N
i
i
i
i
i
i
k
N
k
k
i
n
f
x
l
f
x
lf
x
lf
x
lf
x
lf
x
L
=
+
+
+
+
+
=
=
∑
=
...
...
1
1
0
0
0
.
Каким образом построить базисные полиномы? Определим
( )
(
)(
) (
)(
) (
)
(
)(
) (
)(
) (
)
N
i
i+
i
i-
i
i
i
N
i+
i-
i
x
x
...
x
x
x
x
...
x
x
x
x
x
z
...
x
z
x
z
...
x
z
x
z
z
l
−
−
−
−
−
−
−
−
−
−
=
1
1
1
0
1
1
1
0
, i = 0, 1,..., N.
Легко понять, что
( )
(
)(
) (
)
(
)(
) (
)
( )
(
)(
) (
)
(
)(
) (
)
N
N
N
N
x
x
...
x
x
x
x
x
z
...
x
z
x
z
z
, l
x
x
...
x
x
x
x
x
z
...
x
z
x
z
z
l
−
−
−
−
−
−
=
−
−
−
−
−
−
=
1
2
1
0
1
2
0
1
0
2
0
1
0
2
1
0
,
и т.д.
Функция l
i
(z) является полиномом N-й степени от z, и для
нее выполняются условия «базисности»:
( )
(
)(
) (
)(
) (
)
(
)(
) (
)(
) (
)
N
k
i+
k
i-
k
k
k
N
k
i+
k
i-
k
k
k
k
i
x
x
...
x
x
x
x
...
x
x
x
x
x
x
...
x
x
x
x
...
x
x
x
x
x
l
−
−
−
−
−
−
−
−
−
−
=
1
1
1
0
1
1
1
0
= 0, i
≠
k;
( )
(
)(
) (
)(
) (
)
(
)(
) (
)(
) (
)
1
1
1
1
0
1
1
1
0
=
−
−
−
−
−
−
−
−
−
−
=
N
i
i+
i
i-
i
i
i
N
i
i+
i
i-
i
i
i
i
i
x
x
...
x
x
x
x
...
x
x
x
x
x
x
...
x
x
x
x
...
x
x
x
x
x
l
.
Таким образом, нам удалось решить задачу о построении
интерполирующего полинома N-й степени, и для этого не нужно
решать СЛАУ. Полином Лагранжа можно записать в компакт-
ном виде:
( )
( )
∏
∑
∑
≠
=
=
−
−
=
=
k
i
k
i
k
N
i
i
i
N
i
i
N
)
x
x(
)
x
z(
f
z
lf
z
L
0
0
. Погрешность этой
формулы можно оценить, если исходная функция g(x) имеет
производные до N+1 порядка:
40
]
,
[
),
(
)!
1
(
)
(
)
(
1
0
)1
(
b
a
x
z
N
g
z
r
N
i
i
N
∈
−
+
=
∏
+
=
+
ξ
ξ
.
Из этой формулы следует, что погрешность метода зависит от
свойств функции g(x), а также от расположения узлов интерполя-
ции и точки z. Как показывают расчетные эксперименты, полином
Лагранжа имеет малую погрешность при значениях N < 20. При
бόльших N погрешность начинает расти, что свидетельствует о
том, что метод Лагранжа не сходится (т.е. его погрешность не
убывает с ростом N).
Для вычисления значений полинома нужно выполнить ко-
личество действий, пропорциональное N
2
. В случае небольшого
числа N их можно выполнить на калькуляторе. Если N велико,
для вычислений необходимо использовать ЭВМ. На рис. 3.3
приведен тест программы MathCAD для нахождения в точке
z = 0.35 значения полинома Лагранжа, интерполирующего функ-
цию f = sin(x).
f x()
sin x()
:=
x
0 0.1 0.2 0.3 0.4 0.5 0.6
(
)
T
:=
y
f x()
:=
lagrang z()
s
0
←
p
1
←
p
p
z x
k
−
( )
x
i
x
k
−
( )
⋅
←
i k
≠
if
k 0 6,
∈
for
s
s
y
i
p⋅
+
←
i 0 6,
∈
for
s
:=
lagrang 0.35
(
) =
f 0.35
(
) =
Рис. 3.3. Построение интерполирующего полинома Лагранжа
41
3.5. Метод наименьших квадратов
Во всех вышеизложенных методах условия интерполяции
выполнялись точно. Однако в тех случаях, когда исходные дан-
ные x
i
, f
i
, i = 1,…,N, заданы с некоторой погрешностью
ε
, можно
требовать лишь приближенного выполнения условий интерпо-
ляции: |F(x
i
) – f
i
| <
ε
. Это условие означает, что интерполирую-
щая функция F(x) проходит не точно через заданные точки, а в
некоторой их окрестности, например, как это показано на рис. 3.4.
Приблизим исходные данные глобальным полиномом. Если ре-
шать задачу интерполяции
точно, то полином должен
иметь степень N. При рас-
смотрении полинома Ла-
гранжа мы выяснили, что
полином N-й степени хо-
рошо приближает исход
ную функцию только при
небольших значениях N.
Будем искать полином низкой степени, например,
P
3
(x) = a
1
+ a
2
x + a
3
x
2
+ a
4
x
3
. Если N > 4, то точная задача реше-
ний не имеет: для четырех неизвестных коэффициентов
(a
1
, a
2
, a
3
, a
4
) условия интерполяции дают N > 4 уравнений. Но
теперь точного выполнения условий интерполяции не требуется.
Мы хотим, чтобы полином проходил рядом с заданными точка-
ми. Существует много таких полиномов, каждый из которых
определяется своим набором коэффициентов. Суть метода
наименьших квадратов (МНК) состоит в том, что среди всех
возможных полиномов этого вида выбирается тот, который име-
ет наименьшее среднеквадратичное отклонение в узлах интер-
поляции от заданных значений. Этот многочлен будет самым
близким к заданным точкам из всех возможных многочленов
третьей степени. В i-й точке полином P
3
(x) отклоняется от значе-
ния f
i
на величину (P
3
(x
i
) – f
i
). Суммируем квадраты отклонений
полинома по всем точкам i = 1, 2,…, N, получим функционал квад-
ратов отклонений:
x
y
Рис. 3.4. Приближенное выполнение
условий интерполяции
42
2
3
4
2
3
2
1
1
2
1
3
4
3
2
1
)
(
)
)
(
(
)
,
,
,
(
i
i
i
i
N
i
N
i
i
i
f
x
a
x
a
x
a
a
f
x
P
a
a
a
a
G
−
+
+
+
=
−
=
∑
∑
=
=
.
Найдем минимум этого функционала. Для этого приравня-
ем нулю его частные производные по переменным a
1
, a
2
, a
3
, a
4
.
Используя стандартные правила дифференцирования, получим:
.0
)
(
2
,0
)
(
2
,0
)
(
2
,0
)
(
2
3
4
2
3
2
1
1
3
4
3
4
2
3
2
1
1
2
3
1
3
4
2
3
2
1
2
1
3
4
2
3
2
1
1
=
−
+
+
+
=
=
−
+
+
+
=
=
−
+
+
+
=
=
−
+
+
+
=
∑
∑
∑
∑
=
=
=
=
i
i
i
i
N
i
i
i
i
i
i
N
i
i
N
i
i
i
i
i
i
N
i
i
i
i
i
f
x
a
x
a
x
a
a
x
a
G
f
x
a
x
a
x
a
a
x
a
G
f
x
a
x
a
x
a
a
x
a
G
f
x
a
x
a
x
a
a
a
G
∂
∂
∂
∂
∂
∂
∂
∂
Собирая коэффициенты при неизвестных a
i
, получим
СЛАУ относительно вектора неизвестных (a
1
, a
2
, a
3
, a
4
):
.
,
,
,
3
1
1
4
6
1
3
5
1
2
4
1
1
3
2
1
1
4
5
1
3
4
1
2
3
1
1
2
1
1
4
4
1
3
3
1
2
2
1
1
1
1
4
3
1
3
2
1
2
1
i
N
i
i
N
i
i
N
i
i
N
i
i
N
i
i
i
N
i
i
N
i
i
N
i
i
N
i
i
N
i
i
i
N
i
i
N
i
i
N
i
i
N
i
i
N
i
i
N
i
i
N
i
i
N
i
i
N
i
i
x
f
a
x
a
x
a
x
a
x
x
f
a
x
a
x
a
x
a
x
x
f
a
x
a
x
a
x
a
x
f
a
x
a
x
a
x
a
N
⋅
=
⋅
+
⋅
+
⋅
+
⋅
⋅
=
⋅
+
⋅
+
⋅
+
⋅
⋅
=
⋅
+
⋅
+
⋅
+
⋅
=
⋅
+
⋅
+
⋅
+
⋅
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
∑
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
Полученная система называется нормальной. Для ее реше-
ния используют стандартные методы решения СЛАУ. Как пра-
вило, число неизвестных системы (т.е. число коэффициентов
интерполирующей функции) невелико, поэтому можно приме-
нять точные методы решения СЛАУ, например, метод Крамера
или метод Гаусса.
Метод наименьших квадратов позволяет «приблизить» ис-
ходные данные с помощью линейной комбинации любых элемен-
тарных функций. Часто используются приближения линейной
43
F(x) = a
1
+ a
2
x, тригонометрической F(x) = a
1
sin(x) + a
2
cos(x),
экспоненциальной F(x) = a
1
e
x
+ a
2
e
–x
и т.д. функциями. Некоторые
из них реализованы в MathCAD в виде стандартных функций.
Для
расчета
коэффициентов
линейной
регрессии
y(x) = b + ax можно использовать функцию line(x,y), значением
которой является вектор из двух элементов коэффициентов ли-
нейной регрессии b + ax, а также функции intercept(x,y) и
slope(x,y), возвращающие соответственно коэффициенты b и a
линейной регрессии. Аргументами этих функций являются ис-
ходные данные задачи: x – вектор действительных данных ар-
гумента; y – вектор действительных данных функции того же
размера.
Для
подбора
экспоненциальной
зависимости
вида
2
0
1
)
(
a
e
a
x
E
x
a
+
=
применяется функция expfit(x,y,g), где x и
y – исходные данные, g = (1, 2, 3)
T
. Для подбора степенной за-
висимости
2
0
1
)
(
a
x
a
x
S
a
+
=
используется функция pwrfit(X,Y,g)
с такими же аргументами. Обе функции возвращают вектор, со-
держащий коэффициенты a
0
, a
1
, a
2
.
Применение линейной, экспоненциальной и степенной
функций для аппроксимации данных
x
0 1 2 3 4 5 6
(
)
T
:=
и
y
4.1 3.3 3 4.3 3.6 5.2 5.9
(
)
T
:=
показано на рис. 3.5–3.7. На рис. 3.8
представлены графики линейной
f t()
line x y,
(
)
0
line x y,
(
)
1
t⋅
+
:=
(пунктирная линия), экспоненциальной
f1 t()
exp
0
e
exp
1
t⋅
⋅
exp
2
+
:=
(штриховая линия) и степенной
f2 t()
pw
0
t
pw
1
⋅
pw
2
+
:=
(штрихпунк-
тирная линия) зависимостей в сравнении с исходными данными
(точки). На рис. 3.9 приведены квадраты отклонений аппрокси-
мирующих функций от исходных данных.
x
0 1 2 3 4 5 6
(
)
T
:=
y
4.1 3.3 3 4.3 3.6 5.2 5.9
(
)
T
:=
line x y,
(
) =
f t()
line x y,
(
)
0
line x y,
(
)
1
t⋅
+
:=
Рис. 3.5. Подбор линейной зависимости
44
exp
expfit x y, g,
(
)
:=
g
1
2
3
⎛
⎜
⎜
⎝
⎞
⎟
⎟
⎠
:=
exp
0
=
exp
1
=
exp
2
=
f1 t()
exp
0
e
exp
1
t⋅
⋅
exp
2
+
:=
Рис. 3.6. Подбор экспоненциальной зависимости
pw
pwrfit x y, g,
(
)
:=
exp
0
=
exp
1
=
exp
2
=
f1 t()
exp
0
e
exp
1
t⋅
⋅
exp
2
+
:=
Рис. 3.7. Подбор степенной зависимости
y
f t()
f1 t()
f2 t()
x t,
d
i
d1
i
d2
i
i
Рис. 3.8. График найденных
зависимостей
Рис. 3.9. Графики квадратов
отклонений
45
4. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ
4.1. Численное дифференцирование
Численное дифференцирование, т.е. нахождение значений
производных заданной функции
( )
x
f
y =
в заданных точках x, в
отличие от рассмотренного в п. 4.3 численного интегрирования,
можно считать не столь актуальной проблемой в связи с отсут-
ствием принципиальных трудностей с аналитическим нахожде-
нием производных. Однако имеется ряд важных задач, для ко-
торых численное дифференцирование является единственным
способом нахождения производной. Это, например, поиск про-
изводной таблично заданной функции или дифференцирование
функции в процессе численного решения, когда значения этой
функции известны только в узлах сетки. Кроме того, возможно
сильное усложнение задачи при аналитическом дифференциро-
вании функции, и использование численного подхода упрощает
задачу.
Существует несколько способов для получения формул
численного дифференцирования, которые в конечном итоге мо-
гут привести к одним и тем же формулам. Во-первых, можно
аппроксимировать таблично заданную функцию каким-либо
способом (линейная интерполяция, многочлен Лагранжа,
сплайн-функции и т.д.) и дифференцировать полученную не-
прерывную функцию, приближающую исходную. Эта задача
рассматривалась выше. Во-вторых, для вывода формул числен-
ного дифференцирования можно воспользоваться понятием ко-
нечных разностей. Пусть узлы таблицы x
i
расположены на рав-
ных расстояниях:
ih
x
x
i
+
=
0
, f
i
– соответствующие значения
функции; величину h называют шагом таблицы. Разности f
i+1
– f
i
называют разностями первого порядка. Эти величины обозна-
чают как разность вперед (Δ
+
f
i
), разность назад (Δ
–
f
i+1
): f
i+1
–
– f
i
= Δ
+
f
i
= Δ
–
f
i+1
и центральную разность Δ
±
f
i
= f
i+1
– f
i–1
. Разно-
сти высших порядков образуют при помощи рекуррентных со-
отношений
( )
i
m
i
m
i
m
i
m
f
f
f
f
1
1
1
1
−
+
−
−
Δ
−
Δ
=
Δ
Δ
=
Δ
. Используя эти
формулы, первую производную можно определить разными
способами:
46
( )
h
f
f
x
f
i
i
i
−
=
+
+
1
'
,
(4.1)
( )
h
f
f
x
f
i
i
i
1
'
−
−
−
=
,
(4.2)
( )
h
f
f
x
f
i
i
i
2
1
1
'
−
+
±
−
=
.
(4.3)
Геометрически вычисление производной по трем этим фор-
мулам эквивалентно замене касательной в точке B прямыми BС,
AB и AC соответственно и поиску тангенса угла наклона этих
прямых вместо тангенса угла наклона касательной (рис. 4.1).
x
i–1
f
i–1
A
B
C
Касательная
x
i
x
i+1
f
i–1
f
i+1
Рис. 4.1. Геометрическая интерпретация разностного
дифференцирования
Изучим вопрос о порядке точности (аппроксимации) этих
формул. Разложим f(x) в ряд Тейлора в окрестности точки x
i
:
...
6
2
...,
6
2
3
2
1
3
2
1
+′′′
−′′
+′
−
=
+′′′
+′′
+′
+
=
−
+
i
i
i
i
i
i
i
i
i
i
f
h
f
h
fh
f
f
f
h
f
h
fh
f
f
Подставив эти разложения в (4.1), получаем
( )
...
f
h
f
h
f
h
f
...
f
h
f
h
fh
f
x
f
i
i
i
i
i
i
i
i
i
'
+′′
+′
+′
=
−
+′′
+′
+′
+
=
+
6
2
6
2
2
3
2
Здесь
i
f ′
– первая производная, которую необходимо найти, а
...
6
2
2
+′′′
+′′
i
i
f
h
f
h
– погрешность, с которой вычисляется произ-
водная. Видим, что первый член погрешности имеет порядок h,
значит, при измельчении шага сетки погрешность будет умень-
шаться пропорционально h
1
. Поэтому говорят, что формула
47
имеет первый порядок точности. Нетрудно показать, что фор-
мула (4.2) также имеет первый порядок аппроксимации.
Покажем, что формула центральной разности имеет второй
порядок точности. Подставим в (4.3) разложения для f
i+1
и f
i–1
:
( )
...
f
h
f
h
...
f
h
f
h
fh
f
...
f
h
f
h
fh
f
x
f
i
i
i
i
i
i
i
i
i
i
i
'
+′′
+′
=
+′′
+′
−′
+
−
+′′
+′
+′
+
=
±
6
2
6
2
6
2
2
3
2
3
2
Погрешность вычисления производной пропорциональна h
2
,
значит, формула (4.3) имеет второй порядок аппроксимации.
Основываясь на понятии конечных разностей, можно полу-
чить формулы для аппроксимации производных высших поряд-
ков. Покажем это на примере формулы центральной разности
для аппроксимации второй производной:
( )
2
1
1
2
1
1
2
''
2
h
f
f
f
h
f
f
f
f
h
f
f
x
f
i
i
i
i
i
i
i
i
i
i
−
+
−
+
−
+
±
+
−
=
+
−
−
=
−
=
Δ
Δ
.
Можно доказать, что эта формула имеет второй порядок точно-
сти (доказать самостоятельно).
Наконец, третьим способом получения формул численного
дифференцирования является метод неопределенных коэффи-
циентов. Представим приближенно в точке х = х
0
k-ю производ-
ную таблично заданной функции в виде линейной комбинации
ее значений в узлах
( )
( )
( )
∑
=
≈
N
i
i
i
k
x
f
C
x
f
1
0
,
(4.4)
где C
i
– числовые коэффициенты, которые выберем из условия,
чтобы эта формула была точна для многочлена максимально
высокой степени. Иными словами, потребуем, чтобы для функ-
ции
( )
∑
=
=
m
j
j
j
x
a
x
P
0
приближенное выражение (4.4) было точ-
ным. Для этого необходимо и достаточно, чтобы для любой сте-
пени коэффициенты при a
j
были равны. Поскольку
( )
( )
(
) (
)
k
j
k
j
x
k
j
j
j
x
−
+
−
−
=
1
...
1
, получаем линейную систему
уравнений относительно C
i
:
48
( ) (
)
.
,...,
1,
0
,
1
...
1
0
1
m
j
x
k
j
j
j
x
C
k
j
N
i
j
i
i
=
+
−
−
=
−
=
∑
Если m = N–1, то число уравнений равно числу неизвестных.
Можно показать, что определитель такой системы не обращает-
ся в нуль. Очевидно, что m ≥ k, N ≥ k + 1. При этом точность вы-
числения производной имеет порядок O(h
m+1–k
), хотя при опре-
деленном положении узлов (обычно это симметричное положе-
ние относительно точки x
i
) порядок многочлена, а следователь-
но, и точность можно повысить на единицу.
4.2. Использование стандартных функций MathCAD
для дифференцирования
С помощью MathCAD можно вычислять производные ска-
лярных функций любого количества аргументов: от 0-го до 5-го
порядка включительно. И функции, и аргументы могут быть как
действительными, так и комплексными. Невозможно дифферен-
цирование функций только вблизи точек их сингулярности. Вы-
числительный процессор MathCAD обеспечивает превосходную
точность численного дифференцирования. При этом символь-
ный процессор позволяет вычислять производные громоздких
функций, поскольку, в отличие от всех других операций, сим-
вольное дифференцирование выполняется успешно для подав-
ляющего большинства аналитически заданных функций. В
MathCAD 2001 для ускорения и повышения точности численно-
го дифференцирования функций, заданных аналитически, авто-
матически действует символьный процессор.
Первая производная
Чтобы продифференцировать функцию f(x) в некоторой точке:
1. Определяем точку x, в которой будет вычислена производ-
ная, например, x:=1.
2. Вводим оператор дифференцирования нажатием кнопки
Derivative (Производная) на панели Calculus (Вычисления)
или вводим с клавиатуры вопросительный знак <?>.
3. В появившихся местозаполнителях вводим функцию, зави-
сящую от аргумента x, т.е. f(x), и имя самого аргумента x.
49
Вводим оператор численного <=> или <→> символьного
вывода для получения ответа.
Для численного дифференцирования применяется довольно
сложный алгоритм, вычисляющий производную с точностью:
до 7–8 знака после запятой. Этот алгоритм (метод Риддера) опи-
сан во встроенной справочной системе MathCAD, доступной че-
рез меню Help (Справка). Погрешность дифференцирования не
зависит от констант TOL или CTOL, в противоположность
большинству остальных численных методов, а определяется не-
посредственно алгоритмом. На рис. 4.2 показаны примеры ис-
пользования функций MathCAD для вычисления производных.
x
0.01
:=
x
cos x() ln x()
⋅
d
d
100.041
=
x
cos x() ln x()
⋅
d
d
sin 1. 10
-2
⋅
(
)
−
ln 1. 10
-2
⋅
(
)
⋅
1. 10
2
⋅
cos 1. 10
-2
⋅
(
)
⋅
+
→
Рис. 4.2. Примеры численного и символьного дифференцирования
Производные высших порядков
Чтобы вычислить производную функции f(x) N-го порядка
в точке x, нужно проделать те же самые действия, что и при взя-
тии первой производной, за тем исключением, что вместо опе-
ратора производной необходимо применить оператор N-й про-
изводной (Nth Derivative). Этот оператор вводится с той же па-
нели Calculus (Вычисления) либо с клавиатуры <Ctrl> + <?> и
содержит еще два местозаполнителя, в которые следует помес-
тить число N. «Производная» при N = 0 по определению равна
самой функции, при N = 1 получается обычная первая произ-
водная. Важно перед оператором дифференцирования не забы-
вать присваивать аргументу функции значение, для которого
будет вычисляться производная.
Чтобы вычислить производную порядка выше 5-го, следует
последовательно применить несколько раз оператор N-й произ-
водной подобно тому, как вводятся операторы кратного интег-
рирования. Однако для символьных вычислений этого не требу-
ется – символьный процессор умеет считать производные по-
рядка выше 5-го. Расчет производных высших порядков также
50
производится вычислительным методом Риддера. При повыше-
нии порядка производной на каждую единицу точность падает
примерно на один разряд. На рис. 4.3 показан пример использо-
вания функций MathCAD для вычисления второй производной.
x
0.1
:=
2
x
cos x() x
2
⋅
d
d
2
1.94
=
2
x
cos x() x
2
⋅
d
d
2
1.99 cos .1
( )
⋅
.4 sin .1
( )
⋅
−
→
Рис. 4.3. Пример вычисления второй производной
Частные производные
С помощью символьного и численного процессоров
MathCAD можно вычислять производные функций любого ко-
личества и по разным аргументам, т.е. частные производные.
Чтобы вычислить частную производную, необходимо ввести
оператор производной с панели Calculus (Вычисления) и в со-
ответствующем местозаполнителе напечатать имя переменной,
по которой должно быть осуществлено дифференцирование.
На рис. 4.4 показан пример вычисления частной производ-
ной с помощью стандартных функций MathCAD. В первой
строке определена функция двух переменных, а двух следую-
щих строках символьным образом вычислены ее частные произ-
водные по обеим переменным – x и y – соответственно.
Чтобы определить частную производную численным мето-
дом, необходимо предварительно задать значения всех аргумен-
тов. При этом частные производные высших порядков рассчи-
тываются точно так же, как и обычные производные высших
порядков.
Имеется возможность выбора формы записи частной про-
изводной, причем вид записи не влияет на вычисления, а служит
лишь более привычной формой представления расчетов. Чтобы
изменить вид оператора дифференцирования на представление
частной производной, следует:
1. Вызвать контекстное меню из области оператора дифферен-
цирования нажатием правой кнопки мыши.
51
2. Выбрать в контекстном меню верхний пункт View
Derivative As (Показывать производную как).
3. В появившемся подменю выбрать пункт Partial Derivative
(Частная производная).
Чтобы вернуть вид производной, принятый по умолчанию,
необходимо выбрать в подменю пункт Default (По умолчанию)
либо для представления ее в обычном виде – Derivative (Произ-
водная).
f x y,
(
)
x
2y⋅
cos x() y⋅
+
:=
x
f x y,
(
)
∂
∂
2 x
2y⋅
⋅
y
x
⋅
sin x() y⋅
−
→
y
f x y,
(
)
∂
∂
2 x
2y⋅
⋅
ln x()
⋅
cos x()
+
→
x
1
:=
y
0.1
:=
y
f x y,
(
)
∂
∂
0.54
=
y
f x y,
(
)
∂
∂
cos 1()
→
Рис. 4.4. Пример вычисления частной производной
4.3. Численное интегрирование
Пусть требуется найти значение определенного интеграла
( )
∫
=
b
a
dx
x
f
I
для некоторой заданной на отрезке [a, b] функции
f(x). Хорошо известно, что для функций, допускающих на про-
межутке [a, b] конечное число точек разрыва первого рода, та-
кое значение существует единственно и может быть формально
получено по определению:
( )(
)
1
1
lim
−
=
∞
→
−
∑
=
i
i
n
i
i
n
x
x
f
I
ξ
,
(4.5)
где x
i
– произвольная упорядоченная система точек отрезка [a, b],
ξ
i
– произвольная точка элементарного промежутка [x
i–1
, x
i
].
Из математического анализа хорошо известна формула
Ньютона–Лейбница для нахождения определенного интеграла c
помощью первообразных. Однако первообразную можно найти
не для всех функций, для некоторых элементарных функций
первообразной вообще не существует. Поэтому строятся фор-
52
мулы приближенного интегрирования, которые называются
квадратурными формулами. Простейшие из них выводятся не-
посредственно из определения интеграла, т.е. из представления
(4.5). Зафиксировав там некоторое n ≥ 1, будем иметь
( )(
)
1
1
−
=
−
≈
∑
i
i
n
i
i
x
x
f
I
ξ
.
(4.6)
Это приближенное равенство назовем общей формулой прямо-
угольников. Геометрически площадь криволинейной трапеции
приближенно заменяется площадью ступенчатой фигуры, со-
ставленной из прямоугольников, основаниями которых служат
отрезки [x
i–1
, x
i
], а высотами – ординаты
( )
i
f
ξ
.
Рассмотрим ряд наиболее употребительных квадратурных
формул. Условимся в дальнейшем пользоваться равномерным
разбиением отрезка [a, b] на n частей точками x
i
с шагом
n
a
b
h
−
=
, полагая x
0
= a, x
i
= x
i–1
+ h, x
n
= b. При таком разбиении
формула (4.6) приобретает вид
( )
∈
≈
∑
=
i
n
i
i
f
h
I
ξ
ξ
,
1
[x
i–1
, x
i
].
(4.7)
Теперь дело за фиксированием точек
( )
i
f
ξ
на элементарных
отрезках [x
i–1
, x
i
]. Рассмотрим три случая.
1. Положим
ξ
i
= x
i–1
.
Тогда
из
(4.7)
получаем
( )
( )
∑
∑
−
=
=
−
=
≈
1
0
1
1
n
i
i
n
i
i
x
f
h
x
f
h
I
. Эта формула называется формулой
левых прямоугольников. Ее геометрическая интерпретация при-
ведена на рис. 4.5.
a
b
x
Рис. 4.5. Геометрическая интерпретация формулы
левых прямоугольников
53
2. Пусть
ξ
i
= x
i
. Тогда имеем
( )
∑
=
≈
n
i
i
x
f
h
I
1
.
Это формула правых прямоугольников.
3. Фиксируем
(
)
1
2
1
−
−
=
i
i
i
x
x
ξ
. В результате получаем квад-
ратурную формулу средних прямоугольников:
∑
∑
=
−
=
⎟
⎠
⎞
⎜
⎝
⎛
−
=
⎟
⎠
⎞
⎜
⎝
⎛
+
≈
n
i
i
n
i
i
h
x
f
h
h
x
f
h
I
1
1
0
2
2
.
Покажем, что эта формула имеет второй порядок точности.
Рассмотрим сначала вычисление интеграла на отрезке
[–h/2, h/2],
( )
0
2
2
hf
dx
x
f
h
h
∫
−
≈
,
где f
0
= f(0). Пусть
( )
( )
(
)
2
,
2
1
0
h
F
F
dx
x
f
x
F
x
±
=
=
±
∫
;
( )
( ) (
)
2
2
2
2
h
F
h
F
dx
x
f
h
h
−
−
=
∫
−
.
Разлагая в ряд Тейлора с остаточным членом в форме Лагранжа,
имеем
( )
±
±
±
+
±
=
ξ
''
3
'
2
0
2
1
48
8
2
0
f
h
f
h
f
h
F
,
где
±
ξ
– некоторые точки, такие что
2
2
h
h
<
<
<
−
+
−
ξ
ξ
. То-
гда
( )
( )
2
,
24
''
3
0
2
2
h
f
h
hf
dx
x
f
h
h
<
+
=
∫
−
ξ
ξ
. Для всего интервала [a, b]
формула имеет вид:
( )
(
)
( )
b
a
f
h
a
b
f
h
dx
x
f
n
i
i
b
a
<
<
−
+
=
∑
∫
=
+
ξ
ξ
,
24
''
2
0
2
1
.
54
Таким образом, ошибка численного интегрирования по
формуле средних прямоугольников убывает пропорционально
квадрату шага h, т.е. формула имеет второй порядок точности.
Нетрудно убедиться, что погрешность численного интегрирова-
ния по формулам левых и правых прямоугольников убывает по
линейному закону.
Подстановка в интеграл
( )
∫
b
a
dx
x
f
вместо функции f(x) ее
интерполяционного многочлена Лагранжа той или иной степени
приводит к семейству квадратурных формул, называемых фор-
мулами Ньютона–Котеса. Однако использование в этих фор-
мулах многочленов высоких порядков может быть оправдано
только для достаточно гладких подынтегральных функций. Ча-
ще используются квадратурные правила, получающиеся путем
дробления промежутка интегрирования на большое число мел-
ких частей. Интегрирование на каждой из частей производится с
помощью однотипных простейших формул невысокого порядка.
Приведем два таких правила – трапеций и Симпсона.
Простейшая формула трапеций получается, если на каждом
отрезке [x
i–1
, x
i
] участок кривой (f
i–1
, f
i
) интерполировать линей-
ной зависимостью. Эта формула имеет второй порядок точности
и может быть записана в виде
( )
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+
+
≈
∑
∫
−
=
1
1
0
2
n
i
i
n
b
a
f
f
f
h
dx
x
f
.
Если на каждом отрезке [x
i–1
, x
i
, x
i+1
] участок кривой
(f
i–1
, f
i
, f
i+1
) интерполировать параболой, то придем к формуле
Симпсона, имеющей четвертый порядок точности:
( )
(
)
n
b
a
f
f
f
f
f
h
dx
x
f
+
+
+
+
+
≈
∫
...
4
2
4
3
3
2
1
0
.
4.4. Использование стандартных функций MathCAD
для интегрирования
Интегрирование устроено в MathCAD по принципу «как
пишется, так и вводится». Чтобы вычислить определенный ин-
теграл, следует напечатать его обычную математическую форму
55
в документе. Делается это с помощью панели Calculus (Вычисле-
ния) нажатием кнопки со значком интеграла или вводом с кла-
виатуры сочетания клавиш <Shift> + <7> (или символа «&», что
то же самое). Появится символ интеграла с несколькими местоза-
полнителями, в которые нужно ввести нижний и верхний интер-
валы интегрирования, подынтегральную функцию и переменную
интегрирования. На рис. 4.6–4.8 показаны примеры использова-
ния стандартных функций MathCAD для вычисления интегралов.
0
sin( )
2
x dx
π
=
∫
0
sin( )
2
x dx
π
→
∫
Рис. 4.6. Численное и символьное вычисление
определенного интеграла
α:=2
0
sin( )
4
x dx
π
α
=
∫
x := 1
0
sin( )
42.074
x d
π
α
α
=
∫
Рис. 4.7. Интегрирование функции по разным переменным
g α
( )
0
π
x
α sin x()
⋅
⌠
⎮
⌡
d
:=
i
1 5
..
:=
g i()
2
4
6
8
10
=
Рис. 4.8. Оператор интегрирования в функции пользователя
Чтобы получить результат интегрирования, следует ввести
знак равенства или символьного равенства. В первом случае ин-
тегрирование будет проведено численным методом, во втором –
в случае успеха, будет найдено точное значение интеграла с по-
мощью символьного процессора MathCAD. На рис. 4.6 показаны
оба способа. Конечно, символьное интегрирование возможно
только для небольшого круга подынтегральных функций.
Результат численного интегрирования – это не точное, а
приближенное значение интеграла, определенное с погрешно-
стью, которая зависит от встроенной константы TOL. Чем она
меньше, тем с лучшей точностью будет найден интеграл, но и
больше времени будет затрачено на расчеты. По умолчанию
56
TOL = 0.001. Чтобы ускорить вычисления, можно установить
меньшее значение TOL. Кроме того, пользователь имеет воз-
можность выбирать сам алгоритм численного интегрирования.
Для этого необходимо:
1. Щелкнуть правой кнопкой мыши в любом месте на левой
части вычисляемого интеграла.
2. В появившемся контекстном меню выбрать один из четырех
численных алгоритмов.
При этом возможны четыре численных метода интегрирования:
Romberg – для большинства функций, не содержащих особен-
ностей;
Adaptive – для функций, быстро меняющихся на интервале ин-
тегрирования;
Infinite Limit – для интервалов с бесконечными пределами;
Singular Endpoint – для интегралов с сингулярностью на конце.
Модифицированный алгоритм Ромберга для функций, не опреде-
ленных на одном или обоих концах интервала интегрирования.
Итерационный алгоритм Ромберга применяется, если по-
дынтегральная функция не меняется на интервале интегрирова-
ния слишком быстро и не обращается на нем в бесконечность.
Его основные идеи:
1. Сначала строятся несколько интерполирующих полиномов,
которые заменяют на интервале интегрирования подынте-
гральную функцию f(x). В качестве первой итерации поли-
номы вычисляются по 1, 2 и 4 интервалам. Например, пер-
вый полином, построенный по 1 интервалу, – это прямая
линия, проведенная через две граничные точки интервала
интегрирования, а второй полином – квадратичная пара-
бола и т.д.
2. Интеграл от каждого полинома с известными коэффициен-
тами легко вычисляется аналитически. Таким образом, оп-
ределяется последовательность интегралов от интерполи-
рующих полиномов: I
1
, I
2
, I
4
… Например, по правилу трапе-
ций I
1
= (b – a)
⋅
(f(a) + f(b))/2 и т.д.
57
3. Из-за интерполяции по разному числу точек вычисленные
интегралы I
1
, I
2
, … несколько отличаются друг от друга.
Причем чем больше точек используется для интерполяции,
тем интеграл от интерполяционного полинома ближе к ис-
комому интегралу, стремясь к нему в пределе бесконечного
числа точек. Поэтому определенным образом осуществляет-
ся экстраполяция последовательности I
1
, I
2
, I
4
… до нулевой
ширины элементарного интервала. Результат этой экстрапо-
ляции J принимается за приближение к вычисляемому инте-
гралу.
4. Переход к новой итерации осуществляется с помощью еще
более частого разбиения интервала интегрирования, добав-
ления нового члена последовательности интерполирующих
полиномов и вычисления нового (N-го) приближения Ром-
берга J
N
.
5. Чем больше количество точек интерполяции, тем ближе
очередное приближение Ромберга к вычисляемому интегра-
лу и тем меньше оно отличается от приближения предыду-
щей итерации. Как только разница между двумя последними
итерациями ⏐J
N
– J
N–1
⏐ становится меньше погрешности
TOL или меньше TOL⋅⏐J
N
⏐, итерации прерываются и J
N
по-
является на экране как результат интегрирования.
О расходящихся интегралах
Если интеграл расходится (равен бесконечности), то вычис-
лительный процессор MathCAD может выдать сообщение об
ошибке, выделив при этом оператор интегрирования красным
цветом. Чаще всего ошибка будет иметь тип «Found a number
with a magnitude greater than 10^307» (Найдено число, превы-
шающее значение 10
307
) или «Can’t converge to a solution» (Не
сходится к решению), как, например, при попытке вычислить
интеграл
dx
x
∫
∞
0
1
. Тем не менее символьный процессор справля-
ется с этим интегралом, совершенно правильно находя его бес-
конечное значение.
58
Кратные интегралы
Для вычисления кратного интеграла:
1. Вводится, как обычно, оператор интегрирования.
2. В соответствующих местозаполнителях вводится имя пер-
вой переменной интегрирования и пределы интегрирования
по этой переменной.
3. На месте ввода подынтегральной функции вводится еще
один оператор интегрирования.
4. Точно так же вводится вторая переменная, пределы интег-
рирования и подынтегральная функция (если интеграл дву-
кратный) или следующий оператор интегрирования (если
более чем двукратный) и т.д., пока выражение с многократ-
ным интегралом не будет введено окончательно.
На рис. 4.9, 4.10 приведены примеры символьного и чис-
ленного расчета двукратного интеграла в бесконечных преде-
лах. При этом символьный процессор вычисляет точное значе-
ние интеграла π, а вычислительный определяет его приближен-
но и выдает число 3.142.
∞
−
∞
y
∞
−
∞
x
e
x
2
y
2
+
( )
−
⌠
⎮
⎮
⌡
d
⌠
⎮
⎮
⌡
d
π
→
∞
−
∞
y
∞
−
∞
x
e
x
2
y
2
+
( )
−
⌠
⎮
⎮
⌡
d
⌠
⎮
⎮
⌡
d
3.142
=
Рис. 4.9. Вычисление кратного интеграла
a
b
y
1
−
1
x
x y
3
+
⌠
⎮
⌡
d
⌠
⎮
⌡
d
1
2
b
4
⋅
1
2
a
4
⋅
−
→
a
b
x
1
−
1
y
x y
3
+
⌠
⎮
⌡
d
⌠
⎮
⌡
d
b
2
a
2
−
→
Рис. 4.10. Символьное вычисление кратного интеграла
59
5. ЧИСЛЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
5.1. Численные методы решения задачи Коши
Обыкновенными дифференциальными уравнениями (ОДУ)
можно описывать задачи движения системы взаимодействую-
щих материальных точек, химической кинетики, электрических
цепей, сопротивления материалов (например, статический про-
гиб упругого стержня) и многие другие. Ряд важных задач для
уравнений в частных производных также сводится к задачам для
ОДУ. Так бывает, если многомерная задача допускает разделе-
ние переменных (например, задачи на нахождение собственных
колебаний упругих балок и мембран простейшей формы). Таким
образом, решение ОДУ занимает важное место среди приклад-
ных задач физики, химии и техники.
Рассмотрим ОДУ первого порядка, записанное в общем виде:
( )
y
x
f
dx
dy
,
=
.
(5.1)
Общее решение (5.1) содержит произвольную константу С,
т.е. является однопараметрическим семейством интегральных
кривых
( )
( )
C
dx
y
x
f
x
y
+
=
∫
,
. Для выбора конкретной интеграль-
ной кривой следует определить значение константы С, для чего
достаточно задать начальные данные
y(x
0
) = y
0
.
(5.2)
Несмотря на внешнюю простоту уравнения (5.1), решить
его аналитически, т.е. найти общее решение
( )
C
x
y
y
,
=
с тем,
чтобы затем выделить из него интегральную кривую
( )
x
y
y =
,
проходящую через точку
(
)
0
0
, y
x
, удается лишь для некоторых
специальных типов уравнений. В общем случае решение задачи
можно найти только приближенно.
Приближенное решение задачи (5.1), (5.2) будем искать на
промежутке [x
0
, x
n
]. Построим расчетную сетку
,
0
ih
x
x
i
+
=
i = 1, 2,…,n, с шагом
n
x
x
h
n
0
−
=
. Решение, найденное в узлах
сетки y
i
= y(x
i
), занесем в таблицу
60
x
i
x
0
x
1
…
x
n
y
i
y
0
y
1
…
y
n
На каждом локальном интервале
(
)
1
,
+i
i
x
x
решение задачи
представим в виде
( )
∫
+
=
−
+
1
,
1
i
i
x
x
i
i
dx
y
x
f
y
y
.
(5.3)
Заменим интеграл какой-либо квадратурной формулой чис-
ленного интегрирования. Если воспользоваться простейшей
формулой левых прямоугольников первого порядка точности
( )
(
)
i
i
x
x
y
x
hf
dx
y
x
f
i
i
,
,
1
=
∫
+
,
то получим явную формулу Эйлера:
(
)
i
i
i
i
y
x
hf
y
y
,
1
+
=
+
,
1
...,
,1,
0
−
=
n
i
.
(5.4)
Если в (5.3) использовать формулу правых прямоугольни-
ков, то получим неявный метод Эйлера
(
)
1
1
1
,
+
+
+
+
=
i
i
i
i
y
x
hf
y
y
,
1
...,
,1,
0
−
=
n
i
,
(5.5)
в котором для вычисления неизвестного значения
( )
1
1
+
+
≈
i
i
x
y
y
по известному значению
( )
i
i
x
y
y ≈
требуется решать нелиней-
ное уравнение.
Если для приближенного вычисления интеграла в (5.3) вос-
пользоваться формулой средних прямоугольников
( )
(
)
(
)
2
2
1
/h
x
y,
/h
x
f
h
dx
y,
x
f
i
i
x
x
i
i
+
+
=
∫
+
,
то получим метод предиктор–корректор:
(
)
i
i
i
i
y
x
f
h
y
y
,
2
2
1
+
=
+
,
(
)
2
1
1
,2
/
+
+
+
+
=
i
i
i
i
y
h
x
f
h
y
y
.
(5.6)
Метод (5.6) является частным случаем методов Рунге–Кутты:
(
)
h
y
x
h
y
y
i
i
i
i
,
,
1
ϕ
+
=
+
,
(5.7)
где
61
(
)
( )
∑
=
=
q
n
i
n
i
i
h
k
c
h
y
x
n
1
,
,
ϕ
,
( ) (
)
i
i
i
y
x
f
h
k
,
1
=
,
( )
(
)
h
k
y
h
x
f
h
k
i
i
i
i
1
21
2
,
2
β
α
+
+
=
,
(5.8)
( )
(
)
h
k
h
k
y
h
x
f
h
k
i
i
i
i
i
2
32
1
31
3
3
,
β
β
α
+
+
+
=
,
…………………………………………….
( )
(
)
h
k
h
k
y
h
x
f
h
k
i
q
q
q
i
q
i
q
i
i
q
1
1
,
1
1
...
,
−
−
+
+
+
+
=
β
β
α
.
Здесь
α
n
,
β
nj
, с
n
,
q
n
j
≤
<
<
0
– параметры метода, которые выби-
раются из условия, чтобы метод имел максимально высокий по-
рядок точности p. Можно показать, что p ≤ q. Величины
( )
h
k
i
l
представляют собой правую часть уравнения при различных
значениях аргументов. Их вычисление, занимающее основное
расчетное время, называют этапом метода. Простейшим одно-
этапным методом (q = 1) можно считать метод Эйлера. Наибо-
лее распространенными являются методы Рунге–Кутты второго,
третьего и четвертого порядка (т.е. двух-, трех- и четырехэтап-
ные). При q = 2 получим однопараметрическое семейство двух-
этапных методов Рунге–Кутты второго порядка аппроксимации:
(
)
( )
[
]
.
1
2
,
2
,
<div style="position:absolute;top:73930;left
Информация о работе Рынок: модели, черты