Otvety_ChM_razryv


Чтобы посмотреть этот PDF файл с форматированием и разметкой, скачайте его и откройте на своем компьютере.
1.
Исследование объекта и содержательная постановка задачи


2. Построение математической модели


. Выбор численного метода и разработка вычислительного

алгоритма


4. Разработка программы на компьютере
или выбор пакета
пр
и
кладных програ
мм

5. Проведение вычислений и анализ результатов


1. Математическая модель
,

и
сточники и классификация погрешностей


Процесс решения задачи из физики, техники или экономики с помощью методов
математического моделирования можно представить на следующей схеме:



















Под математи
ческой моделью

понимают систему математических соо
т
ношений
уравнений, неравенств, краевых, началь
ных условий, которым должна удовлетворять
система основных параметров задачи или объекта
. Одно из основ
ных требований
к
математической мо
дели,


соответствие исследуемом
у объекту, т. е.
адек
ватность
.

Другое немалова
ж
ное требование


чтобы модель была
не слишком сложной
,
доступной для мате
матической обработки.

А
налити
ческое решение

наиболее привлекательное,

поскольку позволяют не только
количественно, но и

качественно п
роанализировать исследуемые параметры.
Но
, в
подавляющем большинстве случаев не пред
ставляется во
з
можным

найти анал. реш.
, и


применяются
численные методы
.
Как аналити
ческие, так и численные методы решения
задач подраз
деляются на
точные
и

приближенные
.















Методы решения вычислительных задач

Аналитически
е

Численные

Приближенные

Точные

Алгоритмы

Прямые

Итерационные

О
сновные источники погрешно
стей
:

а

параметры, входящие в описание задачи, заданы неточно; соответству
ю
щую
погрешность называют
неуст
ранимой
;

б

математическая модель описывает изучаемый объ
ект
(
погрешность
матем
ат
и
ческой моде
ли
);

в

численный алгоритм, применяемый для решения математической задачи, зачастую
дает лишь приближен
ное решение
(
погрешность метода
);

г

в процессе вычислений на компьютере промежуточ
ные и конечные р
е
зультаты
округляются
(
вычислитель
н
ая погрешность или погрешность округления
).

П
ервые два вида погрешности, объединяя в один, также называют
неустранимой
погрешностью
.

Обозначив через
I

абсолютную величину погрешности результата, а через
I
Н
,
I
М

и
I
О



абсолютные величины неустранимой погреш
ности, погрешности метода и округ
ления
соответственно


нетрудно получить следующее соо
т
ношение:

I



I
Н

+
I
М

+
I
О
.

(1)

Неравенство 1 дает оценку для погрешности ре
зультата
.

П
олную погрешность результата нельзя сделать ме
ньше, чем наибольшая из
составляющих ее погрешно
стей
.

2. Элементы теории погрешностей

Определение

1.
Приближенным значением

некото
рой величины
а
наз
ы
вается число
а
р
,
которое незначи
тельно отличается от точного значения этой величины.

Определение

2.
Абс
олютной погрешностью

Δ

при
ближенного значения называется
модуль разности меж
ду точным и приближенным значениями этой величины:

Δ

=

|

а

-

а
р
|.

(2)

Определение

3.
Относительной погрешностью

при
ближенной велич
и
ны
а
р

называет
ся отношение
абсолютной погрешности

приближенной вел
и
чины к абсолютной
величине ее точного значения:

δ

=


=
.

(
3
)

Это равенство можно записать в другой форме:

Δ

= |
а
|


δ
.

(4)

На практике, как правило, точное значение величины неизвестно. Поэтому вместо
теоретических понятий абсо
лютной и относительной погрешностей используют прак
-
тические понятия
предельной абсолютной погрешности

и
предельной относительн
ой
погрешности
.

Определение


4.
Под
предельной абсолютной погреш
ностью

прибл
и
женного числа
понимается всякое число
Δ
a
, не меньшее абсолютной погре
ш
ности этого числа:

Δ

=

|

а

-

а
р
|




Δ
а
.
(5)

Неравенство 5 позволяет дл
я точного значения ве
личины получить оценку

а
р

-

Δ
а



а



а
р

+
Δ
а
.
(6)

Часто неравенства 6 записывают в другой форме

а
=
а
р

±
Δ
а

=
а
р

1 ±
Δ
а
). (7)

На практике в качестве предельной абсолютной по
греш
ности выбирают наименьшее
из чисел
Δ
а
,
удовлетворяющих неравенству 5, однако это не всегда возможно.

Определение


5.
Предельной относительной погреш
ностью
δ
а

данного

приближенного числа называется лю
бое число, не меньшее относительной погрешности
этого числа:

δ




δ
а
. (8)

Так как справедливо неравенство

δ 


,

то можно считать, что предельные абсолютная и относи
тельная погрешности связаны
формулой

δ
а

=


или
Δ
а

= |
а
|



δ
а
. (9)

Если абсолютная погрешность Δ
а

значительно меньше точного значения
|а|,
то
относительную погрешность оп
ределяют приближенно как отношение абсолютной по
-
грешности к приближенному значению
:

δ
а




,

Δ
а



|
а
p
|



δ
а
.

(10)

Часто в формуле 10 вместо знака 

ª используют знак точного равенс
т
ва   ª.

Относительную погрешность иногда задают в процентах.
3.
Значащие цифры и правило округления чисел

Определени
е

6.
Значащими цифрами

в записи при
ближенного числа наз
ы
ваются

все
цифры

в его записи, начиная с первой ненулевой слева.

В следующих примерах значащие цифры подчеркнуты.

Определение
7. Первые
п
значащих цифр в записи приближенного числа называются
верными

в узком смысле
, если абсолютная погрешность числа не превосхо
дит половины
единицы разряда, соответствующего
п
-
й

значащей цифре, считая слева направо.

Определение
8. Первые
п
значащих цифр в записи приближенного числа называются
верными
в широком смысле
,
если абсолютная погрешность чи
с
ла не превосхо
дит
единицы разряда, соответствующего
n
-
й

значащей цифре.

Пример

8
.
Пусть
х
 1,1025 ± 0,00009. Верными являются первые чет
ы
ре значащие
цифры, а цифры 5 и  не удовлетворяют определению. В шир
о
ком смысле вер
н
ыми
являются первые пять цифр.

Правило округления чисел

Чтобы округлить число до
п
значащих цифр, отбрасы
вают все цифры, стоящие справа
от
n
-
й значащей цифры, или, если это нужно для сохранения разрядов, заменяют их
нулями. При этом:

1)

если первая отбро
шенная цифра меньше 5, то остав
шиеся десятичные знаки
сохраняют без изменения;

2)

если первая отброшенная цифра больше 5, то к пос
ледней оставшейся цифре
прибавляют единицу;

3)

если первая отброшенная цифра равна 5 и среди ос
тальных отброшенных цифр е
сть
ненулевые, то к после
дней оставшейся цифре прибавляют един
и
цу;

4)

если первая из отброшенных цифр равна 5 и все от
брошенные цифры я
в
ляются
нулями, то последняя остав
шаяся цифра оставляется неизменной, е
с
ли она четная, и
увеличивается на единицу, есл
и нет правило четной ци
ф
ры.

Это правило гарантирует, что сохраненные значащие цифры числа явл
я
ются верными
в узком смысле, т. е. погрешность округлени
я не превосходит половины разря
да,
соответствующего последней оставленной значащей ци
ф
ре. Правило четной

цифры
должно обеспечить ком
пенсацию знаков ошибок.

Пример

10.

Приведем примеры округления до четы
рех значащих цифр:

а ,1415926



3,142;

Δ
р

= |3,142
-

3,1415926| 0,00041 0,0005:

б 1256 410



1256 000;

Δ
р

= |1 256 000
-

1 256 4101
|

500;

в 2
,997925
·

10
8



2,998
·

10
8
;

Δ
р

= |2,998
·

10
8

-

2,997925
·
10
8
| = 0,000075



10
8

0,0005
·

10
8
.

Следующая теорема выявляет связь относительной по
грешности числа с числом
верных десятичных знаков.

Теорема

1. Если положительное приближенное чис
ло имеет
п

верных зн
а
чащих цифр,
то его относительная погрешность
δ

не превосходит велич
и
ны 10
1

-

n
, деленной на первую
значащую цифру
а
н
:

δ



10
1

-

n

/

а
н
. (11)

Формула 11 позволяет вычислить предельную от
носительную погре
ш
ность

δ
a

=

10
1

-

n

/

а
н
. (12)



4.
Погрешности арифметических операций и погрешность вычисления
произвольной функции

Приведем правила вычисления погрешности результа
та различных ари
ф
метических
операций над приближен
ными числами.


Относительно алгебраической суммы
u

=
х
±
у
можно утверждать сл
е
дующее.

Теорема

2. Предельная абсолютная погрешность суммы приближенных чисел равна
сумме предельных абсолютных погрешностей слагаемых, т. е.

Δ
u

=
Δ
x

+
Δ
y
.


(13)

Из формулы 1 следует, что
предельная абсолют
ная погрешность суммы не
может быть меньше предель
ной абсолютной погрешности наименее точного из
сла
гаемых
,
т. е. если в состав суммы входят прибл
и
женные слагаемые с разными
абсолютным
и погрешностями, то сохранять лишние значащие цифры в более точных не

Теорема

. Если все слагаемые в сумме имеют один и тот же знак, то пр
е
дельная
относительная погрешность суммы не превышает наибольшей из предельных относи
-
тельных погрешностей слагаемых
:

δ
u



.
(14)

При вычислении разности двух приближенных чисел
и
=
х
-

у
е
ё

абсолю
т
ная
погрешность, согласно теоре
ме 2, равна сумме абсолютных погрешн
о
стей уменьша
емого
и вычитаемого, т. е.
Δ
u

=
Δ
x

+
Δ
y
,
а предельная от
нос
и
тельная погрешность

δ
u


=
.

(15)

Из формулы 15 следует, что если приближенные значения
х
и
у
близки, то
предельная относительная по
грешность будет очень большой.

Если числа
x

и

y

близки то относител
ьная погрешность будет большой.
В некоторых
случаях удается избежать вычисления разности близких ч
и
сел с помощью
преобразования выра
жения так, чтобы разность была искл
ю
чена
.

Теорема

4. Предельная относительная погрешность произведения
и

=

х


у
приближенны
х чисел, отличных от
н
уля, равна сумме предельных относ
и
тельных
погрешно
стей сомножителей, т. е.

δ
u

=
δ
x

+

δ
y
. (16)

В частности, если
и
=
kx
,
где
k



точное число, имеем
Δ
u

=
|
k
|
Δ
x
,

δ
и

=
δ
х
.

Теорема

5. Предельная относит
ельная погрешность частного равна су
м
ме предельных
относительных по
грешностей делимого и делителя.

Погрешность произвольной функции

Пусть задана произвольная функция
и
=
f
(
x
l
,

x
2
,...,
х
n
),
где
x
l
,

х
2
,


,

х
п



приближенные
величины, а
,
,

,



их известные предельные абс
о
лютные погрешности. Тогда
предельная абсолютная погрешность резул
ь
тата


функции
и


для малых

,
вычисляется

по формуле



(17)

Как видно из формулы 17, для ее применения тре
буется, чтобы фун
к
ция
f
(
x
l
,

x
2
,...,
х
n
)

была дифферен
цируемой по всем переменным.


5.
Аналитический и графический методы локализации корней
уравнения

1
.

Аналитический метод

Теорет
ической основой алгоритма отделения корней служит теорема

К
о
ши

о промежуточных значениях не
прерывной функции.

Теорема

2.1
. Если функция
f
(
x
 непрерывна на отрез
ке [
a
,
b
]

и

f
(
a
) =
A
,

f
(
b
)

=
В
, то для
любой точки
С
, лежа
щей между
А

и
В
, существует точка
c



[
а
,
b
], что

f
(
c
) =
C
.

Следствие
. Если функция
f
(
x
 непрерывна на отрезке [
a
,
b
] и на его ко
н
цах принимает
значения разных зна
ков, то на этом отрезке существует, хотя бы один ко
рень уравнения


f
(
x
) = 0.

Пусть область определения и непрерывности функци
и является коне
ч
ным отрезком
[
a
,
b
]. Разделим отрезок на
n

частей:


a
k

=
а

+
kh
,
k

= 0, 1, ...
п
,
h

= (
b

-

а
)/
п
.

Вычисляя последовательно значения функции в точ
ках
а
0
,
а
1
, ...
,

а
п
, нах
о
дим такие
отрезки [
a
k
,
a
k

+

1
], для которых выполняется условие

f
(
a
k
)

f
(
a
k

+

1
)



0,


(2.1)

т. е.
f
(
a
k
) 0,


f
(
a
k

+

1
  0 или


f
(
a
k
)

� 0,
f
(
a
k

+

1
)



0. Эти отрезки и содержат х
о
тя бы по
одному корню.

Теорема

2.2
. Если непрерывная функция
f
(
x
 монотон
на на отрезке [
а
,
b
] на его
концах принимает значения

разных знаков, то на этом отрезке сущес
т
вует един
ственный
корень уравнения
f
(
x
) = 0.

Если функция
f
(
х
 дифференцируема и ее производная сохраняет знак на отрезке [
а
,
b
], то
f
(
х
 монотонна на этом отрезке.

Если производная
f

'
(
х
)

легко вычисляется и нетру
дно определить ее ко
р
ни, то для
отделения корней уравнения


f
(
х
  0 можно применить следу
ю
щий алгоритм:

1)
найти критические точки, в которых производная
f

'
(
х
 равна нулю или не существует,
и определить ин
тервалы знакопостоянства производной
f

'
(
х
 на
этих интер
валах
функция
f

(
х
 может иметь только по одному корню;

2)
составить таблицу знаков функции
f

(
х
, приравни
вая переменную
х

к кр
и
тическим и
граничным значени
ям, или близким к ним;

3)
определить отрезки, на концах которых функция принимает зна
чения ра
з
ных знаков.

2.
Графический метод

Для отделения корней уравнения естественно приме
нять графический м
е
тод. График
функции
у

=
f

(
х
)

с уче
том свойств функции дает много инфо
р
мации для опре
деления
числа корней уравнения
f

(
х
)

= 0.

С появле
нием мат
ематических п
а
кетов и электронных таблиц ста
ло возможным
вычислять таблицы значений функции с любым шагом и строить графики с высокой
точностью
.

Э
то позволяет уточнять очередной знак в приближенном
зна
чении корня при
помощи следующего алгоритма:

1)
если ф
ункция
f
(
x
 на концах отрезка [
а
,
b
] принимает зн
а
чения разных знаков
,

то делим
отрезок на 10 равных частей

и находим

ту часть, которая содержит ко
рень

таким
способом мы можем уменьшить дл
и
ну отрез
ка
, содержащего корень, в 10 раз;

2)
повторим действия пре
дыдущего пункта для полу
ченного отрезка.

Э
тот процесс можно продолжать до тех пор, пока дли
на

отрезка не станет меньше
заданной погрешности.


6.
Метод половинного деления и метод итераций уточнения корней
уравнения

1.
Метод половинного деления

Метод поло
винного деления 
метод бисекций
,
метод дихотомии
 для р
е
шения
уравнения


f
(
x
  0 заключает
ся в следующем.

Пусть известно, что функция непрерыв
на и принимает на концах отрезка [
а
,
b
] значения
разных знаков, тогда корень содержится в интервале 
а
,
b
).






Метод половинного деления

А
лгоритма метода половинного деления
:

1)
вычислим
х

=

(
a

+
b
)/2;


f
(
x
);

2)
если
f
(
x
  0, переходим к п. 5;

3)
если
f
(
x
)


f
(
a
  0, то
b

=
х
, иначе,
а

=
х
;

4)

если
|
b

-

а
|


ε
, переходим к п.
1;

5)

выводим значение
х
;

6)

конец.

2.
Метод итераций

Метод простых итераций для уравнения
f
(
x
  0 заключается в следу
ю
щем:

1 Исходное уравнение преобразуют к виду, удобному
для

итераций:

x

=
φ
(
х
).



(2.2)

2
 Выбирают начальное приближение
х
0

и вычисляют

последующие прибл
и
жения по
итерационной формуле


x
k

=
φ
(
х
k
-
1
),
k

=1,2, ...


(2.3)

Е
сли существует пред
ел итерационной последовательности
,
,
он является
корнем уравнения


f
(
x
  0, т. е.

f
(
ξ
) =0.







y

=
φ
(
х
)











a

x
0
x
1
x
2

ξ

b




Рис.

2.
4
.

Сходящийся процесс итераций

На рис. 2.
4

показан процесс получения очередного

приближения по мет
о
ду итераций.
Последовательность

приближений сходится к корню
ξ
.

Теоретические основы для применения метода итера
ций дает следующая теорема.

Теорема

2.3
. Пусть

выполняются условия:

1)

корень уравнения
х
=
φ
х
принадлежит отрезку [
а
,
b
];

2)

все значения функции
φ
(
х
 принадлежат отрезку [
а
,
b
],

т. е.
а



φ
(
х
)



b
;

3)

существует такое положительное число
q


1, что производная
φ
'(
x
 во всех точках
отрезка [
а
,
b
] удовлет
во
ряет неравенству |
φ
'(
x
) |


q
.

Тогда:

1)
итерационная последовательность
х
п

=
φ
(
х
п
-
1
)

(
п 
1, 2, , ... сходится при любом
x
0



[
а
,
b
];

2)
предел итерационной последовательности является корнем уравнения


х 
φ
(
x
, т. е. если
x
k

=
ξ
, то
ξ

=
φ
(
ξ
);

3)
справедливо неравенство, характеризующее ско
рость сходимости итер
а
ционной
последовательности

|
ξ
-
x
k
|



(
b
-
a
)

q
k
.

(2.4)


Очевидно что
, эта теорема ставит, довольно, жесткие условия, которые необходимо
проверить перед примене
нием метода и
тераций. Если произво
д
ная фун
к
ции
φ
(
x
 по
модулю больше единицы, то процесс итераций расхо
дится рис. 2.
5
).







Рис. 2.
5
.
Расходящийся процесс итераций



y

=
φ
(
x
)
y

=
x


7.
Метод хорд уточнения корней уравнения

Метод хорд

заключается в замене кривой
у

=
f
(
x
 отрезком

прямой, пр
о
ходящей через
точки 
а
,
f
(
a
)
 и 
b
,
f
(
b
 рис. 2.
6
. Абсцисса точки пересеч
е
ния прямой с осью
ОХ

п
ринимается за очередное приближение.

Чтобы получить рас
четную формулу метода хорд, за
пишем уравнение прямой,
проходящей через точки 
a
,
f
(
a
 и 
b
,
f
(
b
 и, приравнивая
у

к нулю, найдем
х
:





Алгоритм метода

хорд
:

1
)

пусть
k

=

0;

2
)

вычислим следующий номер итерации:
k

=
k

+ 1.

Найд
ем очередное
k
-
e

приближение по формуле:


x
k

=
a
-

f
(
a
)(
b

-

a
)/(
f
(
b
)
-

f
(
a
)).

Вычислим

f
(
x
k
);

3)
если
f
(
x
k
)

 0 корень найден, то переходим к п. 5.

Если

f
(
x
k
)


f
(
b
)



0,
то

b

=
x
k
,
иначе

a

=
x
k
;

4)
если
|
x
k



x
k
-
1
|


ε
, то переходим к п. 2;

5)
выводим значение корня
x
k
;

6)
конец.

Замечание
. Действия третьего пункта аналогичны действи
ям метода пол
о
винного
деления. Однако в методе хорд на каж
дом шаге

может сдвигаться один и тот же конец
отрезка пра
вый или левый, если график функции в о
к
рестност
и корня выпуклый вверх
рис. 2.
6
,
а
 или вогнутый вниз рис. 2.
6
,
б
).

Поэтому в критерии сходимости
используется разность сосед
них приближ
е
ний.





Рис
. 2.
6
.
Метод

хорд





8. Метод Ньютона и

метод секущих
уточнения корней уравнения

1
. Метод Ньютона
(
касательных
)

П
усть найдено приближенное значение корня уравне
ния

f
(
x
)

 0, и об
о
значим его
х
п
.

Расчетная формула
метода

Ньютона
д
ля определения очере
д
ного приближе
ния
x
n
+1

может быть получена двумя способами.

Пер
вый способ выражает геометрический смысл
метода Ньютона

и с
о
стоит в том,
что вместо точки пересечения
гра
фика функции
у
=
f
(
x
)

с осью
О
x

ищем точку
пересе
че
н
ия

с осью
О
x

кас
ательной, проведенной к графику
фун
кции в точке
(
x
n
,

f
(
x
n
)),

как показано на рис. 2.
7
.
уравнение касательной

имеет
вид

у
-

f
(
x
n
)

=
f

'
(
x
n
)(
x

-

x
n
).








Рис. 2.
7
.
Метод Ньютона касательных

В

точке пересечения касательной с осью
О
x

перемен
ная

у
=
0. Приравн
и
вая
у
к нулю,
выразим
х
и получим формулу
метода касательных
:


(2.6)

Приведем достаточные условия сходимости м
етода Ньютона.

Теорема

2.4
. Пусть на отрезке [
а
,
b
]

выполняются ус
ловия:

1)
функция
f
(
x
)

и ее производные

f

'
(
х
)

и
f

''
(
x
)

непре
рывны;

2)
производные
f

'
(
x
)

и
f

''(
x
)

отличны от нуля и сохра
няют определенные постоянные
знаки;

3)


f
(
a
)


f
(
b
)


0 функци
я
f
(
x
)

меняет знак на отрезке.


Тогда существует отрезок [
α
,
β
], содержащий искомый

корень уравнения
f
(
x
)

=
0
, на
котором итерационная пос
ледовательность 2.6 сходится. Если в качестве нулевого
приближения
х
0

выбрать ту грани
ч
ную точку [
α
,
β
], в ко
торой знак функции совпадает со
знаком второй произ
водной,

т.е.

f
(
x
0
)



f
"
(
x
0
)



0, то итерационная последо
вательность сходится моното
н
но рис. 2.
8
).

Замечание
. Отметим, что метод хорд как раз идет с противо
положной стороны, и оба
этих метода могут дру
г друга допол
нять. Возможен и комб
и
нированный
метод хорд
-
касател
ь
ных.



Рис. 2.
8
.
Достаточные условия сходимости метода Ньютона

2
.
Метод секущих


Метод секущих может быть получен из метода Ньюто
на

при замене пр
о
изводной
приближенным выражени
е
м



разностной формулой:


,

,

.



(2.7)

В

формуле 2.7 используются два предыдущих при
бл
ижения

х
п


и

x
n
-
1
.

Поэтому при
заданном начальном
прибли
жении

х
0

необходимо вычислить следующее
приближение

x
1
,
например, методом Ньютона с приближенной
заме
ной производной по формуле

,

Алгоритм метода секущих
:


1)
заданы начальное значение
х
0

и погрешность
ε
. В
ы
чи
с
лим


;

2)
для
п 
1, 2, ... пока выполняется условие
|
x
n



x
n
-
1
|


ε
,
вычисл
я
ем

х
п
+
1
по формуле
(2.7)
.




9.
Метод итераций

решения системы нелинейных уравнений

Система
п
нел
инейных уравнений с
п
неи
з
вестными Имеет вид


f
k
(
x
1
,

х
2
, ... ,

х
п
)
=
0,

1



k



n
.

(2.8)

Приведем систему 2.8 к виду, удобному для итераций:

x
k

=
φ
k
(
x
1
,
x
2
,...,
x
n
, 1 ≤
k


n



(2.11)

Выберем начальное приближение к корню 
, последующие пр
и
ближения
вычислим по формулам

,

s

 0, 1, 2, …

(2.12)

Д
остаточные условия сходимости метода итераций. Об
о
значим точное решение

системы 2.8
.

Определение

2
.3
.
Назовем
ε
-
окрестностью точки

множес
т
во точек
х 
(
х
1
,

х
2
, ... ,

х
п
),

удовлетворяющих условиям

.

Теорема

2.5
. Пусть в некоторой
ε
-
окрестности точно
го

решения

частные производные


существуют и удовл
е
творяют одному из трех
неравенств



(2.13)


где

максимум берется по всем

то
ч
к
а
м

ε
-
окрестности.

Е
сли начальное приб
лижение
х
0

=
(
)

при
надлежит

ε
-
окрестности точного
решения, то метод про
стой

итерации 2.12 сходится к точному р
е
шению.



10.
Метод Ньютона

решения системы нелинейных уравнений

На

практике часто ограничиваются следующим ра
с
суждением
.


Пусть для системы нелинейных уравнений

f
k
(
х
1
,

х
2
, ... ,
х
n
) = 0,

1



k




п
,

в

некоторой
ε
-
окрестности точного решения не равен пулю определитель матрицы
час
т
ных производных мат
рицы
Якоби
):


Тогда существует
начальное приближение
х
0

=(
),
принадлеж
а
щее
ε
-
окрестности точ
ного решения достаточно близкое к точному реш
е
нию, что метод
Ньютона

сходится к точному р
е
шению
.

, (2.14)



k

= 0,1, ...

Для системы двух уравнений

метод
Ньютона

2.14 имеет вид:


(2.15)



k

= 0, 1, ...

11. Нормы векторов и матри
ц

При решении многих прикладных задач весьма полезным являются пон
я
тия нормы
векторов и нормы матриц.

Определение

3
.
1
.
Нормой вектора

называется неотрицательное число, которое
обозначается символом
||
||

и удов
летворяет следующим услов
и
ям:

1)
||
||

 0 при

и
||
||

= 0;

2)
||
||

=
|

c

|

||
||

при любом числовом мн
о
жителе;

3
)

||
||


||
||

+
||
||
. Это соотношение называют неравенством тр
е
угольника.

Норму вектора можно ввести различными способами. На
и
более часто для векторов
n
-
мерного арифметического пр
о
странства


= (
x
1
,
x
2
, ...,
x
n
)
T

используются следующие нормы:

1)
||
||
1
=
,
i

= 1, 2, ...,
n

(3.1)

называется кубической;

2)
||
||
2
=

(3.2)

называется октаэдрической;

3)
||
||
3
=
=

(3.3)

на
зывается сферической, порождена скалярным произведением и определ
я
ет длину
вектора
).

Норма . порождена скалярным произведением которое выражается формулой:
.

Для скалярного произведения векторов справедлив
ы соотнош
е
ния:

=
=
.

Если
A

симметричная матрица, то
=
.

Определение
3
.
2
. Если в пространстве векторов введена норма
||
|
|
, то согласованной
с ней нормой в пространстве матриц называют но
р
му

. (3.4)

Согласованные с нормами векторов .1
-
. нормы матриц определ
я
ются
формулами

,

(3.5)

, (3.6)

. (3.7)

В формуле .7


собственные значения матрицы
A
T
A
, которая явл
я
етс
я
симметр
и
ческой.

Определение

3
.
3
. Две нормы
||
||
α

и
||
||
β

называются
эквивалентными
, если
существуют такие постоянные γ
1

и γ
2
, что для всех векторов

спр
а
ведливы
соотношения:

, и
.

Нормы
||
||
1
,
||
||
2

и
||
||
3

эквивалентны между собой, так как выполнятся нераве
н
ства

||
||
1

||
||
3


|
|
||
2

n

||
||
1
.

Определение

3
.
4
. Будем говорить, что последовательность векторов

сходится к
вектору

по данной норме
||
||
, если выполняется со
о
т
ношение
= 0.

Из эквивалентности норм
||
||
1
,
||
||
2

и
||
||
3

следует, что если последов
а
тельность
векторов сходится по одной из этих норм, то она сходится и по остальн
ым нормам.

П
од нормой матрицы будем понимать норму, согласованную с нормой матрицы.


12.
Метод Гаусса
решения систем линейных алгебраических
уравнений
и его разновидности

Пусть требуется решить систему
п
линейных алгебра
ических уравн
е
ний с
п
неизвестн
ы
ми
:


(3
.
9)

Прямой ход

метода Гаусса преобразует систему .9 к треугольному виду

исключением соответствующих неиз
вестных. Пусть

a
11



0

. Первый шаг заключается в
исклю
чении переменной
x
1

с по
мощью первого уравн
е
ния из ос
тальных уравнений.
Разделим первое уравнение на
a
11
:

;

. (3.10)

Затем от второго уравнения
вычтем

первое уравне
ние, умноженное на
a
21
. В
рез
ультате, на месте второго уравнения получим уравнение, не соде
р
жащее
х
1
.
Чтобы
исключить
х
1

из третьего уравнения,
вычтем из

первое уравнение, умноженное на
a
31
.
Аналогично ис
ключаем
х
1

из четвертого и посл
е
дующих уравнений. Для исключения
х
1

из
i
-
го урав
нения 
i

=

2,

3,

... ,

п
)

приме
ним формулы


;


.

(3.1
1
)


В

результате этих вычислений получим систему вида



(3.12)

На втором шаге исключаем переменную
х
2

с помо
щью
в
торого уравнения из третьего
и последующих уравне
ний.

1.
В общем случае

на шаге
т

для
т

= 1, 2,
...,
n

-

1,

д
елим сначала
т
-
е
уравнение на
:

;

, (3.16)

а
затем исключаем переменную
x
m

c

помощью
m
-
го уравнения из
i
-
го, где

i

=
m

+
1,...,
n
:

;


.

(3.17)

З
десь предполагается, что на каждом шаге выполня
ется условие


0.


В результа
те 
n

-

1)
-
го шага система .9 приобр
е
тает вид


(3.18)

2.

Обратный ход метода Гаусса

вычисляет неизвест
ные
x
i

в обратном порядке. Из
последнего уравнения в .18 нах
о
дим

х
п
=


(
3
.19)

Неизвестные
x
i

определяем по следующим формулам:

x
i

=

i

=
n

-

1,
n

-

2,...,

1.


(3.20)

Метод Гаусса предполагает, что на
т
-
м

шаге выполня
ется условие


0
. Если это
условие не выполняется,

то алгоритм перест
анет работать, так как столкнется с делением
на 0. Кроме этого, в случае выполнения усло
вия


0
, может возникнуть ситуация,
когда ведущий

элемент

бл
и
зок к нулю, что также может привести к неприятностям
в виде

больших п
о
гре
ш
ностей.

Чтобы избежать этих трудностей, применяют
метод Гаусса с выбором главного
элемента
.

Одним из вариан
тов этого метода является
метод Гау
с
са с выбором глав
ного
элемента по столбцам
.

В качестве ведущего эл
е
мен
та на каждом шаге выбираю
т
наибольший по модулю элемент столбца и переставляют соответствующую стро
ку с
другой строкой так, чтобы на
й
денный элемент стал диагональным, затем исключают
соответствующую пере
менную.

Метод Гаусса с выбором главного элемента по стол
бцам

отличается от

алгоритма
(3.16)



.20 только тем, что перед преобразованием .16 ну
ж
но выполнить по
иск
максимального по модулю элемента в
т
-
м

столбце и переставить строки системы
уравнений так, чтобы мак
симальный эл
е
мент стал диагональным элементом матри
цы
к
о
эфф
ициентов.

Приведенный вариант метода Гаусса дает решение и тогда, когда обычный метод
Гаусса пер
е
стает работать из
-
за деления на 0.

В
методе

Гаусса с выбором главного элемента по строкам

на каждом шаге
выбирают наи
б
ольший по мо
дулю элемент строки и перест
авляют столбцы так, чтобы
ведущий элемент находился на диагонали. В этом вар
и
ан
те в уравнениях неизвестные
переменные меняются мес
тами и в алг
о
ритме необходимо думать о том, чтобы запом
-
нить эти пер
е
становки.

1. Решение систем линейных алгебраических у
равнений итерационным
методом и
методом

Зейделя

1
)

Итерационный метод

Запишем систему уравнений .9 в виде

А
х

=
b
,


(3.21)

где
А



матрица коэффициентов;
b



вектор правых ча
стей системы.

Преобразуем .21 к виду, удобн
ому для итераций:


х

 В
х

+
с
.


(3.22)

Тогда метод простой итерации определяется формулой:


x
k
+1

=
B
x
k

+
c
,

k

= 0, 1, ... ,



(3.23)

где начальное приближение
x
0

задано.

В качестве критерия сходимости метода итераций обычно применяют у
с
ловие

||
x
k
+1

-

x
k
||



ε
.

Сформулируем теоремы об условиях сходимости мето
да простых итер
а
ций.

Теорема

3.1

(
достаточное условие сходимос
ти
).

Если ||
В
||  1, то си
с
тема .21
имеет единственное решение, а итерационный метод .2 сх
о
дится к решению со ско
-
ростью геометрической прогрессии.

Теорема

3.2
.

(
необходимое и достаточное условие схо
димости
).

Пусть система
.22 имеет единственное

реше
ние. Итерационный процесс .2 сходится к решению
системы .22 тогда и только тогда, когда все собс
т
вен
ные значения матрицы
В
по
модулю меньше единицы.

Н
а практике для проверки условия сходимости метода итераций более полезна
теорема .1, а теор
ему .2 удается использ
о
вать тогда, когда собствен
ные значения
матрицы
В
известны

Для преобразования системы уравнений к виду, удоб
ному для итераций
, можно
умножить систему .21 на матрицу
D

=
А
-
1

-

ε
, где
ε



произвол
ь
ная матрица. Тог
да
система приме
т вид .22

х

=
В
х

+
с
,


В 
ε
А
,


с

=
D
b


(3.24)

и матрица
В
будет удовлетворять теореме .1, если выб
рать элементы ма
т
рицы
ε

достаточно малыми по модулю. Однако этот прием не выгоден, так как для вычисления
обратной матриц
ы
А
-
1

необходимо выполнить не меньше операций, чем
при

решени
и

самой системы.

При небольшом числе уравнений систему иногда уда
ется привести к виду, удобному
для итераций, с помощью нескольких преобразований.

2
)
Метод Зейделя

Пусть требуется решить систему

уравнений .1:



(3.25)

Выразим из первого уравнения переменную
х
1

из вто
рого


х
2

и т. д.:


Пусть
k
-
e

приближение к решению обозначено как
х
k

=

.
По
д
ставим его в
правую часть полу
ченной системы и выразим следующее пр
и
ближение
x
k
+
1

=


.

В отличие от метода итераций,

метод
Зейделя

использует при
вычислении

уже най
денные компоненты

ве
к
тора


х
k
+1
.

Расчетные формулы метода Зейделя можно записать в виде


(3.26)

Теорема

3.3
(
достаточные условия сходимости
).

Пусть при всех
i

для
коэффициентов системы уравнений .25 выполняются условия

.



(3.27)

Тогда метод Зейделя сходится и выполняется неравен
ство

||
x
k
+1
-

x
k
||
1



q
||
x
k



x
*||
1

(3.28)

где
х
*



точное решение сист
емы .25.

Если матрица
А
удовлетворяет условию .27,
гово
рят, что
А


матр
и
ца с
диагональным преобладанием.

Теорема

3
.
4
(
достаточные условия сходимости
).

Пусть матрица
А

си
с
темы .1


вещественная, симмет
ричная и положительно определенная. Тогда ме
тод Зей
деля
сходится.


14.
Погрешность решения и обусловленность системы уравнений

Рассмотрим влияние погрешности правой части и
с
войств матрицы си
с
темы
линейных уравнений на по
гр
ешность решения. Пусть правая часть системы задана
пр
иближенно, с погрешнос
тью
η
:

А
х

=
b
1
,

b
1

=
b

+
η
.

Пусть
х
1

решение неточно заданной системы
А
х

=
b
1
,


х



решение то
ч
ной системы
А
х

=
b
.

Обозначим по
гр
ешность решения через
r

=
х
1

-

х
.

Тогда мо
ж
но запи
ть
А
х
1

=
b
1

в виде
А
(
х

+
r
)

=
b

+
η
,

и
А
r

=
η
.

Определение

3
.
5
.
Мерой обу
словленности системы

на
зывается число

. (3.29)

Мера обусловленности системы равна верхней грани
от
ношения относ
и
тельной
погрешности решения к отно
си
тельной погрешности правой части. Из формулы .29

сл
еду
ет неравенство


(3.30)

Если мера обусловленности системы принимает боль
шое значение, это значит, что
небольшая погрешность правой части может привести к бол
ь
шой погрешности
решения, т. е.
полученное приближенное решение ока
жется непригодным.

Учитывая, что
r

 А
-
1
η
,

можно получить формулу вы
числения меры
обусловленности системы:


(3.31)


Определение

3
.
6
.
Мерой обусловленности матрицы

А
называется чи
с
ло


.

(3.32)

Для вычисления меры обусловленности матрицы мож
но с помощью .1 получить
формулу


(3.33)

Учитывая .0, можно записать



(3.34)

Неравенство .4 связывает относительные погреш
ности правой части и решения
системы через свойства матрицы системы.

Определение

3
.
7
. Системы уравнений и матрицы назы
ваются
плохо об
у
словленными
,

если их меры обусловл
ен
ности принимают большие знач
е
ния, и
хорошо обуслов
-
ленными
,

если меры обусловленности принимают м
а
лые значения.

Понятно, что при решении хорошо обусловленных систем малые п
о
грешности
правой части приводят к ма
лым погрешностям решения, а пл
о
хо об
у
словле
нные си
-
стемы уже нельзя решать обычными методами.

15. Вычисление определителя и обратной матрицы с помощью метода
Гаусса

Вычисление определителя матрицы является класси
ческим примером з
а
дач, для
решения которых важно найти эффективные алгоритмы.

При неп
осредственном раскрытии определителя квад
ратной матрицы

n
-
го порядка
нужно найти сумму
n
 слага
емых, каждое из которых равно произведению
n

элементов
матрицы, взятых по одному с каждого столбца и каждой строки, т. е. нужно выполнить
порядка
n
!
n

операций.


Если матрица приведена к диагональному или треу
гольному виду, то ее
определитель равен произведению диагональных элементов.

Для преобразования матрицы к треугольному виду можно применить м
е
тод Гаусса,
который

потребует порядка 2
n
3
/ операций.

Если вним
ательно проанализировать метод Гаусса, то можно увидеть, что он
фактически позволяет одновре
менно с приведением матрицы коэффиц
и
ентов к треу
-
гольному виду, вычислить определитель этой матрицы. Дейс
т
вительно, определитель
матрицы коэффициентов системы урав
нений .18 равен произведению диагональ
ных
элементов, т. е. единице. С другой ст
о
роны,
если к любой строке матрицы прибавить
другую строку, умно
женную на число, определитель не изменится. А если строку
матрицы разделить на число, отличное от нуля, то о
пределитель матрицы
разделится на то же чи
с
ло.

Отсю
да следует, что в результате преобразований к виду
.18, опред
е
литель матрицы коэффициентов системы уравне
ний .9 изменился на
множитель, который равен про
изведению ведущих элементов, т. е. мы пол
у
ч
аем форму
-
лу для вычисления определителя


Метод Гаусса также эффективен и для вычисления обратной матрицы.
Вычисление
обратной матрицы мож
но рассматривать как частный случай решения совокупно
-
сти систем линейных уравнений с одной и
той же матри
цей коэффициентов и
разными правыми частями
. В этом случае преобраз
о
вания, ко
торые применяются к
столбцу
правых частей системы, нужно пр
и
менить к несколь
ким столбцам правых
частей.


16. Собственные числа и собственные векторы
матрицы,

и их
в
ычисление

Определение

3
.
8
. Собственным числом или собствен
ным значением квадратной
матрицы
А
называется чис
ло
λ
, такое, что система уравнений

А
х

=
λ
х


(3.35)

имеет
ненулевое решение

х
.
Это решение называется соб
ственным вектором матриц
ы
А
,

соответствующим соб
ственному значению
λ
.

Собственный вектор определяется с точностью до по
стоянного множит
е
ля
,
если
х

удовлетворяет .5, то и
с
х
также является решением .5.

Преобразуем систему .5 к виду
(
А
-

λ
Е
)

х

=
0
, где
Е


единичная
ма
т
рица. Так как
система линейных од
нородных уравнений имеет ненул
е
вые решения лишь тогда, когда
определитель матрицы равен нулю, получим уравнение для определения собственных
значений


(
A



λ
E
) = 0, (3.36)

кото
рое называется
характеристическим
,

или
вековым
,

уравнением
.

Е
сли матрица диагональная или треугольная, определ
и
тель равен произведению
диаго
нальных элементов и поэтому собственные зн
а
чения рав
ны диагональным
элементам. Нетрудно вычислить соб
ственные зна
чения для трехдиагональной матрицы, а
также для почти треугольной матрицы.

Для диагональной матрицы собственному значению
λ
i

=
a
ii

отвечает ед
и
ничный
собственный вектор
х
i

= (0, ... , 1, ...
,0)
T
,
у которого
i
-
я

компонента равна 1, а остальные
компоненты 0
.

Теорема

3
.
5
. Собственные значения симметричной матрицы с действ
и
тельными
элементами действительны, а собственные векторы, соответс
т
вующие различным соб
-
ственным значениям, взаимно ортогональны.

Теорема

3
.
6
.

Если
λ
min

и
λ
m
ах



наименьшее и наиболь
шее соб
ственные значения
действительной симметричной матрицы
А,
то для любого вектора
х

справедливо
неравен
ство

λ
min
(
x
,

х
)



(
А
х
,

х
)



λ
max
(
х
,

х
).


(3.37)

Определение

3
.
9
.

Действительная симметричная мат
рица
А
называется положительно
опред
еленной, если для любого вектора
х



0

выполняется у
с
ловие

(
А
х
,

х
)


0.


(3.38)

Теорема

3
.
7
.

Действительная симметричная матрица
А
является полож
и
тельно
определенной тогда и только тогда, когда все ее собственные знач
е
ния положитель
ны.

Теорема

3
.
8
(
критерий Сильвестра
).

Для того чтобы действительная симметричная
матрица
А 

(
а
ij
)

была по
ложительно определенной, необх
о
димо и достаточно, что
бы все
главные диагональные миноры ее определ
и
теля были положительны:



(3.39)

Теорема

3
.
9
(
теорема Перрона
).

Если все эл
е
менты квадратной матрицы
положительны, то ее наибольшее по модулю собственное значение полож
и
тельно и не
являет
ся кратным, а соответствующий собственный ве
к
тор име
ет положительные
ко
ординаты.

Рассмотрим
итерационный метод

определения наи
большего по модулю
собственного значения и соответ
ствующего собственного вектора матрицы
А
,

который
запишем в виде следующего алгоритма.


Алгоритм определения наибольшего по модулю соб
ственного

знач
е
ния и
соответствующего собственного вектора матрицы с полож
и
тельными
элементами
:

1)
зададим начальное приближение
х
0

к собственному вектору;
k

= 0;

2)
в
ы
числим следующие приближения
х
k
+1

по формулам

k
+1

=
A
x
k
,

,
x
k
+1

=
,



k

=

k
+
1;


(3.40)

 если
|
λ
k
+1

-

λ
k
|


ε
,

переходим к п. 2, иначе


к п. 4;

4)

конец.

Критерием для остановки итераций является условие
|
λ
k
+1

-

λ
k
|
ε



з
а
данная
погрешность.

С помощью формулы 
.40 можно вычислить снача
ла
k
-
ю

степень ма
т
рицы
А

и
умножить ее на вектор
х
0

см. пример

10, а для
λ
k

можно брать о
т
ношение ненулевой

координаты вектора

х
k
+1

к соответствующей координ
а
т
е

вектора
х
k
,
которая также не
должна быть равной нулю. Так как зар
анее н
е
известно, какие координаты соб
ственного
вектора не равны нулю, можно брать отноше
ние сумм координат, если эта сумма не
равна нулю.


Метод скалярных произведений

Рассмотрим метод скалярных произведений для опре
деления наибольш
е
го
собственного знач
ения и соответ
ствующего собственного вектора дейс
т
вительной
матри
цы
А
.

Теорема

3
.
10
.

Транспонированная матрица
А
T

имеет те же собственные значения, что
и матрица
А.
Пусть
λ
i

и
λ
k



различные собственные знач
е
ния матрицы
А
и
транспонированной матрицы
А
T
)
,

х
i



собственный вектор матрицы
А,
отвечающий
собственному значению
λ
i

а
y
k



собственный ве
к
тор матрицы
А
T
,
отвечающий
собственному значению
λ
k
.

Тогда векторы
х
i

и
у
k



орто
гональны.

Пусть требуется вычислить наибольшее собственное значение и соотве
т
ст
вующий
собственный вектор дей
ствительной матрицы
А
.

В
методе ск
а
лярных произведе
ний

вместе с матрицей
А
используется транспонирован
ная матрица
А
Т
.

Алгоритм метода скалярных произведений
:


1 з
ададим начальные приближения:
х
0



к

собствен
ному
в
ектор
у матр
и
цы
А
и
у
0

=
х
0



к собственному век
тору транспонированной матрицы
А
T
;


k

= 0;


2)
вычислим 
k

+
l
)
-
e

приближение к наибольшему собственному знач
е
нию
λ

по
формулам

x
k

+
1

=

А
х
k
,



y
k

+1

=
A
T
y
k
,

λ
k

=
,

k

=
k

+1;




(3.41)


3)
если
|
λ
k
+
1

-

λ
k
|



ε
, переходим к п. 2, иначе


к п. 4;


4 к
онец.

Вычисление всех собственных значений

положительно определенной
симметричной

матрицы

Приведем алгоритм для вычисления нескольких пер
вых или
всех собс
т
венных
значений и соответствующих собственных векторов положительно определенной сим
-
метричной матрицы.

Пусть уже вычислены первые
m

собственных значений
λ
1
,
λ
2
,
...,
λ
m

и

т

соответствующих собственных векторов

x
1
,
x
2
, ...,
x
m
.


Алгоритм вы
числения очередного
(
т

+

1)
-
го

собственного значения и

с
о
ответствующего

собственного вектора
:


0)
выберем начальное приближение
;
k

=

0;

1)
вычислим
k
-
e

прибл
и
жение к собственному значе
нию
λ
m
+1
:

;

(3.42)

2 находим вектор
из уравнения

; (3.43)

3)


если
m

 0, ортогонализируем вектор

к
первым

т
собственным векторам:


-

; (3.44)

4)


нормируем полученный вектор:

; (3.45)

5)

k

=
k
+
1.

Процесс 1

5 повторяется до тех пор, пока не
будет выполнено условие сходимости
итераций



(3.46)

где
ε



заданная погрешность.

При вычислении первого собственного значения и со
ответствующего вектора п. 
пропускается.

Этим алгоритмом можно в
ычислить все собственные значения и собс
т
венные
векторы.


17.
Интерполяция и
интерполяционные формулы Ньютона

Под задачей
интерполирования
,

или
интерполяции

принято понимать следующую
:

Пусть известны значения
y
i

функции
f
(
x
)

в точках
х
i

i

= 0, 1,
...,
п
,

называ
е
мых

узлами
интерполяции
.

Требу
ется найти такую функцию
F
(
x
)

(
инте
р
по
лирующая фун
кция
),

значения которой в узлах интерпол
я
ции совпада
ют со знач
е
ниями
f
(
x
):

F
(
x
i
) =

y
i
,

i

=
0,

l
, ... ,
n
.


(4.1)

Если зависимость и
н
терполирующей фун
кции
F
(
x
)

от параметров лине
й
на, то говорят
о
линей
ной интерполяции
,

в противном случае


о
нелинейной интерполяции
.

Формула
у
=
F
(
x
)

называется
интерполяционной

и используется для в
ы
числения
значения функции
f
(
x
)

в точке
х
,

не совпадающей с узлами инте
р
п
оляции.
Эта операция
называется
интерполированием интерполяци
ей функции
.

При этом, если точка
х
не
принадл
е
жит от
резку
[
a
,

b
],

a

=
min
(
x
0
,
х
1
,

... ,

х
п
),

b

=

m
ах
x
0
,
х
1
,

... ,

х
п
)
,

то говорят об
экстраполировании фун
к
ции
.

Часто в качестве интерполирующ
ей функции исполь
зуется многочлен
т
а
кого

порядка
Р
п
(
х
)

(
интерполяционный многочлен
),

удовлетворяющий у
с
ловию


y
i

=
P
n
(
x
i
)
,

i

= 0,
l
,...,
п.


(4.2)

Так как многочлен
n
-
го порядка определяется своими коэффициентами, то число
параметров равн
о
n

1 и ус
ловия 4.2 представляют систему лине
й
ных уравн
е
ний
.


(4.3)

Определителем системы 4. является о
п
ределитель Вандермонда


(4.4)

который о
тличен от нуля, если среди узлов интерполяции нет совпада
ю
щих. Отсюда
следует, что задача построения интерполяционного многочлена ра
з
решима и существует
единственный интерполяционный многочлен, удовл
е
тво
ряющий условиям 4.2. Как
увидим ниже, формы запи
с
и интерполяцио
н
ного многочлена могут быть разными.
Отметим, что задача отыскания
пр
о
извольного много
члена
,

удовлетворяющего
условиям 4.2, некорректна. Е
с
ли степень многочлена меньше
п
,

то такого многочле
на
не существует, но если больше
п
,

то таких мног
очле
нов бесконечное множ
е
ство.


Интерполяционные формулы Ньют
о
на

Пусть узлы интерполяции распределены на отрезке равномерно
х
i

=
а

+
ih
,
i

= 0, 1, ... ,
п
.

Обозначим шаг из
менения переменной
х

ч
е
рез

х

=
h
.

Введем сначала понятия
конечной разности

и
обоб
щенной степени
,

к
о
торые
используются для записи ин
терполяционной формулы Ньют
о
на.

Первой конечной разностью

функции
у
=
f
(
x
)

называ
ется выраж
е
ние


у
=

f
(
x

+


х
)



f
(
x
).


(4.5)

Конечная разность второго порядка

опр
е
деляется формулой


2
у
=

(

y
)

=

[
f
(
х

+

х
)
-

f
(
х
)]
=
f
(
x

+
2

х
)
-

2
f
(
х

+

х
) +
f
(
x
). (4.6)

Используя формулу бинома Ньютона, можно вывести формулу
конечной разности п
-
го п
о
рядка
:


n
у
=

(

n
-
1
y
)

=
f
(
х

+

n

х
)

-

f
[
х

+

(
n

-

1)

x
]

+
f
[
x

+

(
n



2
)

х
]

-

...+
(
-
1)
n

f
(
x
)
,
где
. (4.7)

Справедливо утвержде
ние:
конечная разность

(
n
+
1)
-
го

порядка мног
о
члена п
-
го
п
о
рядка равна нулю
.

Пусть функция
f
(
x
)

задана своими значения
ми
y
i

=

f
(
x
i
)

в равноотстоящих точках
x
i

=

х
0

+
ih
,

i

= 0, 1,
...
,

п
.

Таблица конечных разностей

функции
f
(
х
 запис
ы
ва
ется в одной из двух форм:
горизонтальной
или

диаго
нальной таблицы разностей
.
Так как к
о
нечная разность
использует 2 значения, столбец

y
i

содержит
п
значений,

2
y
i



на одно значение
меньше и т. д. Если задано
п

+
1 значений фун
к
ции, то таблица конечных разностей со
-
держит
п
столбцов, причем последний столбец содержит только одно значение. Общий
вид
горизонтал
ь
ной таб
лицы
,
конечных разнос
тей

прив
е
ден в табл. 1.

Таблица

1


x
i


y
i


y
i


2
y
i


3
y
i

...


n
-
1
y
i


n
y
i


x
0


y
0


y
0


2
y
0


3
y
0

...


n
-
1
y
0


n
y
0


x
1


y
1


y
1


2
y
1


3
y
1

...


n
-
1
y
1


...

...

...

...

...

...



x
n
-
2


y
n
-
2


y
n
-
2


2
y
n
-
2





x
n
-
1


y
n
-
1


y
n
-
1






x
n


y
n








Обобщенной степенью п
-
го порядка

числа
х
называется произвед
е
ние

x
[
n
]

=
x
(
x
-

h
)(
x
-

2
h
)(
x
-

3
h
)
...[
x


(
n


1)
h
], (4.8)

где
h

=
const
,
а для
n

= 0
полагают
x
(0)

= 1.



Найдем конечные разнос
ти для обобщенной

степени:


х
[
n
]

=
(
х 
h
)
(
n
)

-

x
(
n
)

=
(
x

+
h
)
x
(
x

-

h
)(
x

-

2
h
)...[
x

-

(
n

-

2)
h
]

-

x
(
x

-

h
)(
x

-

2
h
)...[
x

-

(
n

-

1)
h
]

=

x
(
x

-

h
)(
x

-

2
h
)(
x

-

3
h
)...[
x

-

(
n

-

2)
h
]{(

x

+
h
)

-

[
x

-

(
n

-

1)
h
]}

=
nh

x
[
n

-

1]




Мы получили следующу
ю формулу:


х
[
n
]

=

nh

x
[
n

-

1]


(4.9)

Пользуясь методом математической
индукции, мож
но
вывести и общую формул
у



k
х
[
n
]

=
n
(
n

-

1)...[
n

-


(
k

-

1)]

h

k
x
[
n

-

k
]
.

(4.10)

Выведем
первую интерполяционную формулу Ньютона
.
Най
дем мног
о
член

Р
п
(
х
),
удовлетворяющий условиям
Р
п
(
х
i
)

=
y
i

для


i

= 0, 1,...,
n
.

Будем и
с
кать

многочлен

Р
п
(
х
)

в
сл
е
дующем виде:


P
п
(
х
)

=
а
0

+
a
1
(
х



х
0
) +
a
2
(
x



x
0
)(
x



x
1
)

+ ...
+
a
n
(
x



x
0
)(
x



x
1
)...(
x

-

x
n
) =

=
а
0

+
a
1
(
х



х
0
)
[1]

+
a
2
(
x



x
0
)
[2]

+ ..
. +
a
n
(
x



x
0
)
[
n
]
.

И
з условия
Р
п
(
х
0
)

=
y
0

следует
a
0

=
y
0
.

Найдем первую конечную разность

многочлена
Р
п
(
х
):



Р
n
(
х
) =
a
1
h

+
2
h
а
2
(
x

-

x
0
)
[
1
]

+ ... +
n
a
n
(
x

-

x
0
)
[
n
-
1
]
.

Отсюда при
х

=
х
0

полу
чим

Р
n
(
х
0
)

=
a
1
h

=

y
0
, т.е.
a
1

=

y
0
/
h
.
Аналогично можно


по
луч
ить общую формулу

коэффициентов многочлена Ньютона:

a
k

=
. (4.11)

Теперь можно записат
ь

первую

интерполяционную формулу Нь
ю
тона
:

P
п
(
х
)

=
y
0

+
(
х

-

х
0
)
[1]

+
(
x

-

x
0
)
[2]

+ ... +
(
x

-

x
0
)
[
n
]

(4.12)

С помощью замены переменной
q

= (
х
-

x
0
)
/
h

первую интерполяционную формулу
Ньютона можно предста
вить в виде

P
п
(
х
)

=
y
0

+
q

+

+ ... +
.

(4.13)

Аналогичными рассуждениями выводится
вторая интерполяционная формула
Ньют
о
на
:

P
п
(
х
)

=

y
n

+

(
х

-

х
n
)

+
(
x

-

x
n
)
(
x

-

x
n
-
1
)

+..
.
+
(
x

-

x
n
)
(
x

-

x
n
-
1
)
...
(
x

-

x
1
)
,(4.14)


которая с помощью замены
q

=
(
х
-

x
n
)/
h

приводится к в
и
ду

P
п
(
х
)

=

y
n

+

q

+

+ ... +
.

(4.1
5
)

Первая интерполяционная формула Ньютона приме
няется для интерпол
и
рования
вблизи начала таблицы около
х
0
, а вторая


для интерполиров
а
ния вблизи кон
ца
таблицы около
х
n
).


18.
Интерполяционная формула Лагранжа
.
Остаточный

член
интерполяционной формулы Лагранжа

Интерполяционные формулы Ньютона 4.12, 4.15 пригодны лишь для
равноотстоящих узлов интерпол
и
рования. Рассмотрим интерполяционную формулу
Лагран
жа, которая применима при любом расположении узлов. Пусть дана систе
ма
точек

х
0
,

х
1,

...
,

х
n
,
принадлежащих некоторому отрезку [
а
,
b
], и известны
соответствующие значения функции
y
i

=
f
(
x
i
).

Найдем мн
о
гочлен
L
n
(
x
),

удовлетворяющий
услов
и
ям
L
n
(
х
i
) =
y
i
.

При построении многочлена Лагранжа используются многочлены
n
-
й ст
е
пени

р
i
х,

принимающие значение 1 в точке
x
i

и нулевые значения в остал
ь
ных узлах интер
поляции
x
j
,

j



i
.

Так как
x
j
,

при


j



i

являются

корнями мн
о
гочлена


р
i
(
х
)
,
то справедливо
разложение

р
i
(
х
)


на мно
жители

р
i
(
х
)

=
C
i
(
х

-

х
0
)(
х

-

х
1
)...(
х

-

х
i
-
1
)
(
х

-

х
i
+
1
)

...
(
х

-

x
n
),

i

= 0, 1, ... ,
n
.

Из условия
p
i
(
x
i
)

=
1 находим значение ко
н
станты
C
i
:


C
i
=
1/
(
х
i

-

х
0
)(
х
i

-

х
1
)...(
х
i

-

х
i
-
1
)
(
х
i

-

х
i
+
1
)

...
(
х
i

-

x
n
)

и получаем выражение для многочленов
p
i
(
x
i
):

р
i
(
х
)

=
.
(4.1
6)

Интерполяционный многочлен Л
а
гранжа

имеет вид

L
n
(
x
) =

=



(4.17)

Формулу Лагранжа можно преобразовать так, чтобы упростить вычисл
е
ния. Вынесем
за знак суммы произве
дение

П
n
+1
(
x
)
=

(
х

-

х
0
)(
х



х
1
) ... (
х

-

х
n
),

которое является многочленом степени
п
1. Получим

L
n
(
x
) =
П
n
+1
(
x
)
.

Теперь формулу Лагранжа можно записать в в
и
де

L
n
(
x
) =
П
n
+1
(
x
)
, (4.18)

где
D
i

=
(
х
i

-

х
0
)
(
х
i

-

х
1
)...(
х
i

-

х
i
-
1
)
(
x

-

x
i
)
(
х
i

-

х
i
+
1
)

...
(
х
i

-

x
n
)
.

Обратите внимание на то, что в произведении
D
i

из
n

+

1 сомножителя вместо скобки
(
х
i

-

x
i
)

присутствует множитель
(
х


x
i
).

Остаточный член интерполяционной формулы Лагранжа

Для анализа погрешности
интерп
оляции использует
ся остаточный

член
и
н
терполяционного многочлена Лаг
ранжа

f
(
x
)
-

L
n
(
x
) =
=

(4.19)

где


принадлежит отрезку
[
a
,

b
],
который определяется числами

а
=
min
(
x
0
,
х
1
,
...,
х
п
,
х
),

b

=
m
ах
x
0
,
х
1
,
...,
х
п
,
х
).


(4.20)

Формула 4.19 выводится с помощью функции

φх
) =
f
(
x
)
-

L
n
(
x
)



K
П
n
+1
(
x
),

где параметр
К
в
ы
бирается из условия
φ
(
х
) = 0:

K
=
(
f
(
x
)
-

L
n
(
x
)
)/
П
n
+1
(
x
).

19.
Равномерное приближение функции. Многочлены Чебышева

Оп
ределение

4
.
1
. Многочлен
Р
п
(
х
)

равномерно прибли
жает

на отре
з
ке
[
а
,

b
]
фун
кцию

f
(
x
)
с точностью до
ε
, если выполняется нер
а
венство


|

f
(
x
)

-

P
n
(
x
)|



ε
.

(4.21)

Т
еорем
а Вейерштрасса:

Теорема

4
.
1
.

Если функция
f
(
x
)

непрерывна на отрез
ке
[
а
,

b
],

то для люб
о
го
ε

� 0
найдется многочлен
Р
n
(
х
 достаточно высокой степени
n
, который равномерно при
-
ближает на отрезке
[
а
,

b
]

функцию
f
(
x
)

с точностью до
ε
, т. е. выполняется 4.21.

Определение

4
.
2
. Многочлен
Р
п
(
х
)

называется
много
членом наилуч
шего
приближения

для функции
f
(
x
)

на отрезке
[
а
,

b
],

если для любого многочл
е
на
Q
n
(
x
)
степени
n

выполняется нераве
н
ство

|

f
(
x
)

-

P
n
(
x
)|



|

f
(
x
)

-

Q
n
(
x
)

|.



(4.22)

Многочлены Чебышева

оп
ределяются рекуррентными формулами

Т
0
(
х
) = 1,

Т
1
(
х
) =
х
,

T
n
+
1
(
x
) =
2
хТ
п
(
х
)

-

Т
п
-
1
(
х
)

при
n




1.

(4.23)

Выпишем многочлены Чебышева
Т
п
(
х
 для
п
= 2, 3, 4, 5:

Т
2
(
х
)

= 2
хТ
1
(
х
)
-

T
0

= 2
х
2

-

1,

Т
3
(
х
) =
2
хТ
2
(
х
)

-

T
1

=
4
х
3

-

3
х
,

Т
4
(
х
) = 8
х
4

-

8
х
2

+ 1,

Т
5
(
х
) = 16
х
5

-

20
х
3

+ 5
х
.

Старший коэффициент многочлена
Т
n
(
х
, т. е. коэффи
циент при
х
п
,

равен
2
n
-
1
.

Справедливо представление многочленов Чебышева через тригонометр
и
ческие
функции

T
n
(
x
)

=
cos
(
n

arccos

х
, при
п


0.


(4.24)

Поэтому

|
Т
n
(
х
)|



1

для

|

х

|



1.


(4.25)

Корни многочлена Чебышева принадлежат отрезку

[
-
1; 1]:

x
k

=
,
k

=
0
, 1, ... ,
n

-

1.



(4.26)

Многочлены Чебышева
Т
п
(
х
 с четными индекса
ми являются четными функциями, а с
нечетн
ы
ми


нечет
ными функциями.

Многочлены Чебышева
Т
n
(
х
 обладают следующим замечательным сво
й
ством. Если
умножить
Т
п
(
х
)

на 2
1
-
n
, получится многочлен,
наименее укл
о
няющийся от нуля

на
отрезке [
-
1; 1].

Теорема

4
.
2
. Если

Р
п
(
х
)



произвольный многочлен степени
п
со ста
р
шим
коэффициентом 1, справе
д
ливо не
равенство

|

P
n
(
x
)

|



|

2
1
-
n

T
n
(
x
)

| =

= 2
1
-
n

,

(4.27)

где

=
2
1
-
n

T
n
(
x
).

В

некоторых учебниках по численным методам многочленом Чебышева называется
мног
о
член

Т
n
(
х
) =
2
1
-
n

T
n
(
x
)
,

наименее уклоняющийся от нуля

на отрезке [
-
1; 1].

Поясним смысл теоремы 4.
2

на примерах.

Если
п
 2, то теорема 4.
2

утверждает, что наибольшее значение любой квадратной
функции вида
х
2

+
рх 
q

на отрезке [
-
1; 1] не меньше
2
1
-
2

= 0,5.

Наибольшее значение любого многочлена вида
х
3

+

рх
2

+

qx

+

r

на отре
з
ке [
-
1; 1] не
меньше 2
1
-
3

= 0,25.

Среди всех квадратных функций вида

х
2

 рх 
q

наи
менее уклоняющимся от нуля на
отрезке [
-
1; 1] много
членом является функция

2
1
-
2

T
2
(
x
)

=

х
2

-

0,5.

Среди всех многочленов вида
х
3

+

рх
2

+

qx

+

r


наиме
нее укл
о
няющимся от нуля на
отрезке [
-
1; 1] является многочлен

2
1
-
3
Т
3
(
х
)

=

(4
х
3

-

3
х
)/4

=
х
3

-

0,75
х
.

Для произвольного отрезка
[
а
,

b
]

многочлен со стар
шим коэффициентом 1,
наименее укл
о
няющийся от

нуля, получается из

заменой

.

Этот многочлен имеет вид

= (
b

-

a
)
n

2
1
-
2
n

T
n

Корнями многочлен
а


являются точки

x
k

=
,
k

=
0
, 1, ... ,
n

-

1.

(4.29)

Многочлены Чебышева используются для минимиза
ции погрешности
и
н
терполяционной формулы за счет

оптимального выбора узлов интерполяции.
В
формуле остаточного члена 4.19 у многочлена
П
n
+1
(
х
)

старший коэфф
и
циент равен 1.
Для минимизации погрешности интерполяционной формулы для функции
f
(
x
)
на отрезке
[
а
,

b
]

нужно взять в качестве узлов интерпол
я
ции
точки

x
k

=
,
k

=
0
, 1, ... ,
n
.

(4.
30
)

являющиеся корнями многочлена
. Тогда по
грешность интерпол
я
ции
оценивае
т
ся неравенством

|
f
(
x
)
-

L
n
(
x
)
|


.

(4.31)




20. Аппроксимация. Метод наименьших квадратов в случае линейной
аппроксимирующей функции

Наиболее распространенным методом аппроксимации

построение кривой которая
была бы максимально близка к экспериментальным точкам
)

экспериментал
ь
ных

данных
является
метод наимень
ших квадратов
.

Пусть заданы значения функции
y
i
,

соответствующие значениям
x
i
,
i

= 1, 2, ... ,
п.
Предположим, что
аппрок
симирующая

функция

g
(
x
)

зависит от
т
параметров

g

=
g
(
x
,

a
1
,
a
2
,
... ,

a
m
),

т



п
.

При
точечной

квадра
тич
ной а
п
проксимации

параметры
а
1
,

а
2
,

... ,
а
т

аппрокси
мирующей функции опред
е
ляются из условия минимума суммы квадратов
отклонений значений а
п
проксимирую
щей фун
к
ции от заданных значений функции:

,

(4.51)

Вид функции
g
(
x
,
a
1
,

а
2
,

... ,
а
т
)

определяется особен
ностями решаемой з
а
дачи
.

Если аппроксимирующая функция линейно зависит от параметров, то м
е
тод
наименьших квадратов приводит задачу ее определения к системе л
и
нейных уравнений.

Наиболее часто применяются аппроксимации прямой
линией
(
линейная регрессия
),

полиномом
(
полиномиаль
ная регрессия
),

линейной комбинац
и
ей линейно независи
-
мых функций
.

Рассмотрим применение метода наименьших квадратов для определения параметров
функции
g
(
x
,
a
1
,

а
2
, ... ,
а
т
)

при аппроксимации
линейной ко
мб
и
нацией линейно не
-
зависимых функций
:

g
(
x
,
a
1,

а
2
, ... ,
а
т
)

= a
1
φ
1
(
x
) +

a
2
φ
2
(
x
)

+
...+
a
m
φ
m
(
x
).


(4.52)

Тогда условие 4.51 имеет вид

S
(
a
1
,...,
a
m
) =
.


(4.53)

Минимум квадратичной функции существует; необхо
димым услов
и
ем его
существования явл
я
ется равенство нулю частных производных

,

k

=
1, ... ,
m
.

(4.54)

Из 4.54 получим систему уравнений метода наимень
ших квадратов
н
е
обходимого
ус
ловия существования минимума:



(4.55)

При
линейной аппроксимирующей функции

g
(
x
,

а
0
,

а
1
) =

а
0

 а
1
х


(4.56)

система 4.55 имеет вид




(4.57)

Аналогично можно по
лучить систему уравнений для определения пар
а
метров
полиномиальной регрессии:

g
(
x
,

а
0
, ... ,
а
т
) =

а
0

+

а
1
х
+ ... +

а
т
х
т
.


(4.58)

В этом случае получим систему из
т

+
1 уравнения с
т
1 неи
з
вестным


(4.59
)

Если
зависимость

аппроксимирующей функции от параметров нел
и
нейна
,
то для
определения параметров приходится решать нелинейные зад
а
чи. Ин
о
гда удается с
помощью преобразований вместо функции с нелиней
ной зависимостью от параметров
рассматривать функцию

с линейной зав
и
симостью от параметров.




21.
Графическое дифференцирование

Пусть известен график функции
у 
f
(
х
 на отрезке
[
а
,
b
].

Можно построить график
производной функции, вспомнив ее геометрический смысл. Воспользуемся тем фактом,
что производ
ная функции в точке
х
равна тан
генсу угла наклона к оси абсцисс
касательной к ее графи
ку в этой точке.

Если
х

=

х
0
,

найдем
у
0

=

f
(
x
0
)

с

помощью графика и затем проведем касательную
АВ
к
графику функции в точке 
х
0
,

y
0
 рис. 5.1. Проведем прямую, паралл
ельную касательной
АВ,
через точку 
-
1,

0 и найдем точку
у
1

ее пересечения с осью ординат. Тогда значение
у
1

равно тан
генсу угла наклона касательной к оси абсцисс, т. е. про
изводной функции


f
(
x
)

в точке
х
0
:

у
1
=

=
tg
α

=
f


(
x
0
)
,
и точка
М
0

(
х
0
,
у
1
)
принадлежит графику
производной.







Рис. 5.1

Чтобы построить график производной, необход
имо разбить отрезок
[
а
,

b
]

на
несколько частей точками
х
i
, затем для каждой точки графически построить значение
производной и соединить полученные точки плавной кри
вой с помощью лекал.

На рис. 5.2 показано построение пяти точек
М
1
,

М
2
,

... ,
М
5

и графика
производной.

Алгоритм построения графика производной
:

1.


Строим касательную к графику функции
у
=
f
(
x
)

в точке
(
х
1
,

f
(
x
1
));

из точки 
-
1, 0)
параллельно касатель
ной в точке
(
х
1
,

f
(
x
1
))

проведем прямую до пересечения с осью
ординат; эта точка пересечения дае
т значение про
изводной
f


(
х
1
).

Строим точку
М
1
(
х
1
,
f


(
х
1
)
)
.

2.


Аналогично построим остальные точки
М
2
,

М
3
,
М
4

и
М
5
.

3.


Соединяем точки
М
1
,

М
2
,

М
3
,

М
4
,

М
5

плавной кривой.

















M
4



Рис. 5.2

Полученная кривая является графиком производной.

Точность графического способа определения производ
ной невысока. Мы приводим
описание этого способа толь
ко в учебных целях.

Замечание
. Если в алгоритме построения графика производ
ной вместо точки 
-
1,

0)
взять точку
(
-
l
,

0, где
l

 0, то график будет построен в другом масштабе по оси ординат.





22. Разностные формулы для обыкновенных производных

Разностные формулы для приближенного вычисления

производной подсказаны самим
опр
еделением производной.

Пусть значения функции в точках
x
i

обозначены через
y
i
:

y
i

=
f
(
x
i
),

x
i

=
a
+
i
h
,

i

=
0,
1, ... ,
n
;
h

=


Мы рассматриваем случай равномерного распределения точек на отрезке
[
a
,

b
]. Для
приближенного вычи
сления производных в точках
x
i

можно использовать следующие
разностные формулы
,
или
разностные производные
.

. (5.1)

. (5.
2
)


. (5.3)


Так как предел отношения 5.1 при
h



0 равен пра
вой производной в точке
х
i
, то это
отношение иногда на
зывают
правой разностной производной

в точке
x
i
.

По
аналогичной причине отношение 5.2 называют
левой разностной производной

в точке
x
i
.

Отношение 5. на
зывают
центральной разностной производной

в точке
x
i
.

Оценим погрешность разностных формул 5.1

5., предполагая, что функция
f
(
x
)
разлагается в ряд Тейло
ра в окрестности точки
x
i
:

f
(
x
)

=
f
(
x
i
)

+
. (5.4)

Полагая в 5.4
х
=
x
i

+
h

или
х  х
i

-

h
,
получим

y
i
+
1

=
f
(
x
i

+
h
) =
y
i

+
. (5.5)

y
i
-
1

=
f
(
x
i

-

h
) =
y
i

-

.
(5.
6
)


Учитывая 5.5 и 5.6, име
ем


,


(5.
7)


,


(5.8)


.




(5.9)


Из последних соотношений следует, что разностная формула 5. имеет погрешность
на порядок меньшую, чем разностные формулы 5.1 и 5.2.

Производные высших порядков можно приближенно вычислять по формулам,
получе
нным с помощью после
довательного применения разностных соотношений 5.1

(5.3).

Разностная формула для второй производной
(
разно
стная производная второго
порядка
)

имеет вид


. (5.10)


Н
епосредственной подстановкой разложений 5.5 и 5.6 в формулу 5.10
можно
получить зависимость между второй производной функции и
разностной формулой
для производной второго порядка
:

. (5.11)




2. Ра
зностные формулы для частных производных

Разностные формулы для частных производных ана
логичны формулам 5.1

(5.3),
(5.10).

Пусть функция двух переменных
и
=
f
(
x
, у
)

определе
на в прямоугольной области
а
1



х


b
1
,

а
2



у



b
2
.


Определение
5.1
.
Назовем
сеточной областью

множе
ство точек
(
x
i
,

y
j
),

где

x
i

= a
1

+
ih
1
,

i
= 0, 1, ... ,
n
1
;


h
1

=
(
b
1

-

a
1
)
/
n
1
,

y
j

=
a
2
+
j
h
2
,


j =
0,

l, ... ,

n
2
;


h
2

=
(
b
2

-

a
2
)
/n
2
.

На рис. 5.
4

изображена сеточная область для
п
1

= 5,
п
2

=
5, которая состоит из 6
точек.




Введем обозначение
u
ij

=
f
(
x
i
,

y
j
).

Тогда для частных производных первого и второго порядка по переменной
х
можно
записать разностные формулы 5.1

5. и 5.10:

.

(5.1
2
)

.

(5.
13
)

.

(5.
14
)



(5.
15)


Аналогичные разностные формулы можно записать и для частных производных
первого и второго порядка по перемен
ной
у
:

.

(5.1
6
)

.

(5.
17
)

.

(5.
1
8
)



(5.
19)



Запишем разностные формулы для смешанных

производных:


(5.20)


(5.21)


(5.22)




24. Вычисление производных с помощью интерполяционных формул
в
случае равномерного распределения узлов

Пусть заданы зн
ачения функции
y
i

=

f
(
x
i
)

в узлах
х
i
,

равномерно распределенных на
отрезке
[
а
,

b
]:

x
i

=
а

+
ih, i
= 0, 1,
...
,

n
;

h =
(
b
-

a
)
/n
.

Построим формулы приближенного дифференцирова
ния с помощью
первой
интерполяционной формулы Нью
тона

4.1, которую после некот
орых упрощений
мож
но записать в виде

Р
n
(
х
)

=
y
0
+
q

y
0

+
, (5.23)

где
q

= (
x



x
0
)
/
h
.

Мы предполагаем, что узлы интерполяции распреде
лены на отрезке
[
а
,

b
]

равномерно
с шагом
h
.

Выберем в качестве
х
0

узел, ближайший к точке
х.
Диффе
ренцируя 5.2 по
переменной
х
и учитывая, что

,

получим приближенные формулы для
вычисления про
изводных:

, (5.24)

.

(5.2
5
)

Чем больше слагаемых в формулах 5.24 и 5.25 и
с
пользуется, тем выше точность
вычисления производ
ных. Эти формулы дают хорошие приближения для то
чек, близких
к значению
х
0

левому концу отрезка
[
а
,
b
]).

Если
х  х
0
,

то
q

=
0 и для вычисления производных в узлах
x
i

получим формулы

,

(5.2
6
)

.

(5.2
7
)

Очевидно, что
вторая интерполяционная формула Ньютона

позволяет вывести
формулы для вычисления производных в точках, близких к правому краю
х
n

отрез
ка
[
а
,

b
]:


Р
n
(
х
)

=
y
n
+
q

y
n
-
1

+
, (5.2
8
)

,

(5.2
9
)

.

(5.
30
)

При
x

=
x
n

получим

,

(5.
31
)

.

(5.
3
2
)

Выведем формулы для вычисления производной с по
мощью интерполяционного
мног
очлена Лагранжа:

L
n
(
x
) =


,


y
(
x
)
=

L
n
(
x
)
+



В случае равномерного распределения узлов с помо
щью обозначений

q

=
,

q
[
n
+1]

=
q
(
q

-
1)...(
q

-

n
)

получим

П
n
+1
(
x
) =

(
х
-

х
0
)(
х



х
1
) ...

(
х

-

х
п
)

=

h
n

+

1
q
(
q

-
1)...(
q
-
n
)

=
h
n
+
1

q
[
n
+1]
,

П
׳
n
+1
(
x
i
) =
(
х
i

-

х
0
)(
х
i



х
1
) ...

(
х
i

-

х
п
)

=
h
n
i
(
i
-
1)...1 (
-
1)...(
-
n
+
i
)

=
(
-
1)
n
-
1
h
n
i
!(
n
-
i
)!.

Тогда полином Лагранжа запишем в виде

L
n
(
x
)

=
.

Найдем производную

y
׳
(
x
) =
L
׳
n
(
x
)

=

. (5.33)

С помощью 5. получим формулы для вычисления производной в точках
x
i

при
различных значениях
п
:

1)


п

 2. Используются три точки, и производная в

этих точках выражается формулами



(5.34)


2)


n

  четыре точки.



(5.35)




3)


n

= 4 (
пять

точек
).




(5.36)





25.
Вычисление производных с помощью интерполяционных формул
в
случае неравномерного распределения уз
лов

Пусть известны значения
y
i

функции в узлах
х
i
:

x
0


x
1


x
2

...
x
n
-
1

x
n
.

В этом случае для вычисления производных использу
ют интерполяционный
многочлен. Предположим, что точка
х
расположена ближе к начальному узлу
х
0
,

тогда
мы можем применить
первый

интерполяционный много
член Ньютона

y
(
x
) =
y
(
x
0
) +(
x

-

x
0
)
y
(
x
0
,
x
1
) + (
x

-

x
0
)(
x

-

x
1
)
y
(
x
0
,
x
1
,
x
2
) + (
x

-

x
0
)(
x

-

x
1
)(
x

-

x
2
)
y
(
x
0
,
x
1
,
x
2
,
x
3
)
+ ...

Здесь использованы обозначения для
разделенных раз
ностей
:

,

,

,...

Обозначим через

разность
(
х
-

x
i
)

и запишем много
член Ньютона в виде

y
(
x
) =
y
(
x
0
) +

ξ
0
y
(
x
0
,
x
1
) +
ξ
0
ξ
1
y
(
x
0
,
x
1
,
x
2
) +
ξ
0
ξ
1
ξ
2
y
(
x
0
,
x
1
,
x
2
,
x
3
) + ...
(5.37)

Теперь можем вывести формулы для производн
ых:

у'
(
х
)

 ух
0
,

х
1
)

+
(
ξ
0

+

ξ
1
)
ух
0
,

х
1
,

х
2
)
+ (
ξ
0
ξ
1
+
ξ
1
ξ
2

+

ξ
0
ξ
2
)

y
(
x
0
,
x
1
,
x
2
,
x
3
) + ...

у"
(
х
)

=
2

y
(
x
0
,
x
1
,
x
2
)

+
2(
ξ
0

+
ξ
1
+

ξ
2
)

y
(
x
0
,
x
1
,
x
2
,
x
3
) + ...

y
(
k
)
(
x
)

=

k
![
y
(
x
0
,

x
1
,

...
,

x
k
)

+

(
ξ
0

+

ξ
1

+
...
 ξ
k
)
y
(
x
0
,

x
1
,
...
,

x
k
+1
)

+

+

(5
.
38)

Оставляя в 5.8 несколько слагаемых, получим фор
мулы для приближенного
вычисления производных. При этом порядок погрешности по отношению к шагу разбие
-
ния
h

равен числу оставленных членов, или разности меж
ду числом узлов ин
терполяции
и порядком производной.


Приведем некоторые простые формулы 
h

=
max

h
i
):

1. Первая производная по двум точкам:

y
'
(
x
)

=
y
(
x
0
,
x
1
)
+
O
(
h
)

=



(5.39)


Первая производная по трем точкам:

,

, (5.40)


2. Вторая производная по трем точкам



(5.41)


Вторая производная по четырем точкам:

[

]
.


(5.42)

Третья производная по четырем точкам:

.




26.
Практическая оценка погрешности вычисления производных. Метод
Рунге
-
Ромберга

Рассмо
трим метод Рунге для повышения порядка точ
ности формул для вычисления
производных на примере формулы 5.11:


(5.44)

Предполагая, что функция
f
(
x
)

четырежды дифферен
цируема, погрешность
разностной фор
мулы для второй производной можно представить в виде


(5.45)

Вычислим по формуле 5.44 вторую производную в той же точке
х
i
, но используя
сетку с шагом
k
·
h
,

k



целое, например с
k

= 2:


(
5.46)


Погрешность формулы 5.46 имеет вид


(5.47)

Вычитая 5.45 из 5.47, получим

-


(5.48)

Отсюда для погрешности получим
первую формулу Рунге
:


(5.49)

Теперь вторую производную в точке
x
i

можно вычис
лить по уточненной формуле
(
второй формуле Рунге
):


(5.50)




27. Численное интегрирование. Квадратурные формулы Ньютона
-
Котеса

Формулы для приближенного вычисления интегралов часто, называют
квадр
а
турными.

Пусть требуется вычислить определе
н
ный интеграл


b







y
(
x
)
dx

=
F
(
b
)

-

F
(
a
)



(6.
2
)


а

При выводе квадратурных формул для приближенно
го вычисления опр
е
деленного
интеграла вспомним его геометрический смысл


инте
грал равен площади кри
-
волинейной трапеции, ограниченной граф
и
ком подынтег
ральной функции, осью
Ох

и
отрезками прямых
х

=
a

и
х

=
b
.

Разобьем отрезок [
а
,

b
] на
п

частей точк
а
ми
х
i
:


x
i

=
a

+
ih
,
i

= 0, 1,...,
n
;
x
0

=
a
,
x
n

=
b
,

(6.3)

Обозначим через
y
i

значения функции в точках
x
i
. За
меним подынтеграл
ь
ную
функцию интерполяционным многочленом Л
а
гранжа 4.17:

y
(
x
)

=
L
n
(
x
) =

=




(
6
.
4
)

Тогда получим приближенную фо
рмулу для вычисле
ния интегр
а
ла
:







,




(6.5)

где
A
i



числовые коэффициенты, которые не зависят от поды
н
тегральной функции и их
значения для заданного
п

всегда
можно опр
е
де
лить.

Выведем формулы для вычисления коэффициентов
A
i
. Введем обознач
е
ния

,
q
[
n

+

1]

=
q
(
q

-

l
)...(
q

-

n
. Тогда многочлен Лагранжа можно зап
и
сать в виде

L
n
(
x
) =

,

Заменяя под знаком интеграла в 6.5 функцию
у
(
х
 многочл
е
ном
L
n
(
x
),

получим


=
,



=
.


Отсюда следуют формулы для вычисления коэфф
и
ци
ентов
A
i
:

A
i

=

=

=
·

= (
b



a
)
H
i
,

где

i

= 0, 1,...,
n
.

Коэффициенты
H
i

называются коэффиц
и
ент
ами Коте
ca

и

вычисляются по фо
р
мулам

H
i

=
,
i

= 0, 1,...,
n
. (6.6)




28.
Формула прямоугольников и

формула трапеций

1
)

Формула прямоугольников

Формулы прямоугольников получаются заменой по
дынтегрально
й фун
к
ции
постоянным значением. В каче
стве такого знач
е
ния выбирают значение функции в од
ной
из точек отрезка
[
а
,

b
]

или на левом конце отрезка, или на правом конце, или в сер
е
дине
отрезка рис. 6.1:





y
(
a
)(
b



a
)
,



(6.
7
)





y
(
b
)(
b



a
)
,



(6.
8
)





y
(
)(
b



a
)
,



(6.
9
)






Если формулы 6.
7
)

(6.
9
 применить к каждой части
[
x
i
,
x
i

+

1
]

частичного
отрезка
[
а
,

b
],

получим
общие формулы прямоу
гольников
.

Фактически, о
п
ределенный интеграл
в этом
случае
прибли
женно заменяется интеграл
ь
ной су
м
мой:





,


(6.10)





,


(6.11)





,


(6.12)

Геометрически это означает, что площадь криволиней
ной трапеции пр
и
ближенно
з
а
меня
ется площадью сту
пенчатой фигуры. В частности, рис. 6.2
иллюстрир
у
ет формулу
(6.10).







Рис. 6.2.

Формулу 6.10 называют
формулой левых прямоуголь
ников
, (6.11)



фо
р
мулой правых
прямоугольников
, а

(6.12)



формулой средних прямоугольн
и
ков
.

Формулы прямоугольников практически не использу
ются из
-
за большой погрешности
порядка
O
(
h
 у форму
лы средних прямоугольников более выс
о
кий порядок
O
(
h
2
)).

2
)

Формула трапеци
й

Положим в формул
е

(6.6)

п

 1 и вычи
с
лим значения
H
i
:

H
i

=
,
i
= 0, 1;

H
0

=

=

=

=
;


H
1

=

=
,
откуда

A
0

=
A
1

=
.

Таким образом,
замен
яя

подынтегральную функцию многочле
ном Лагра
н
жа первой
ст
е
пени получ
аем

формулу тра
пеций
:




h
,
h

=
b



a
.



(6.
13
)

Геометрический смысл формулы трапеций 6.1 зак
лючается в том, что кривая
у

=

у
(
х
)
заменяется отрезком прямой, проходящей ч
е
рез точки 
х
0
,

у
0
 и 
х
1
,

у
1
, или, в др
угих
обозначен
и
ях, 
а
,
у
(
а
 и
(
b
,
у
(
b
 рис. 6..

Понятно
, что формулы трапеций и средних прямоу
гольников являются точными для
линейной фун
к
ции.

Если обобщить 6.1 для равномерного разбиения от
резка на
п

частей, приходим к
общей формуле трапеций

р
ис. 6.4:





=
,
.



(6.
1
4
)


y

x

a

b

0

Рис. 6.






Рис. 6.
4



Погрешность

формулы трапеций 6.1 есть в
е
личина поря
д
ка
O
(
h
3
. В этом можно
убедиться, используя фор
мулу погрешности интерполяционного мн
о
гоч
лен
а Лаг
ранжа. А
для общей формулы тр
а
пеций 6.14 погреш
ность есть величина порядка
O
(
h
2
, так как при
сумми
ровании погрешности накаплив
а
ются.




29.
Формула Симпсона

Применяя интерполяционную формулу Лагранжа при
п

 2, получим зн
а
чения
коэфф
и
циентов:

H
i

=
,
i
= 0, 1, 2.

H
0

=

=

,

H
1

=
-

=
,

H
2

=

=
,

h

=
,
A
i

= (
b



a
)
H
i

= 2
hH
i
,

A
0

=
,
A
1

=
,

A
2

=
.




. (6.15)

Геометрический смысл формулы Сим
п
сона 6.15 зак
лючается в том, что кривая
у

=
у
(
х
 заменяется частью парабол
ы, проходящей через три точки 
х
0
,
у
0
), (
х
1
,
у
1
 и 
х
2
,
у
2
)
рис. 6.5, а.

Формула Симпсона точна не только для полинома вто
рой степени, но и для полинома
третьей степени в силу симметрии, показа
н
ной на рис. 6.5, б.

Остаточный член формулы 6.15 имеет п
орядок
O
(
h
5
).

Общая формула

Симпсона
строится для четного
п

= 2
m
; при этом формула 6.15

применяе
т
ся для каждо
го отре
з
ков
[
х
0
,
х
2
], [
х
2
,
х
4
], ... , [
х
п

-

2
,
x
n
]:





(
у
0

 4у
1

+ 2
у
2

+ 4
у
3

+
...+
2
у
2
m
-
2

+ 4
y
2
т
-
1

+
у
2
т
) =

=

(6.16)

x

x

y

y

0

0

x
0

x
1

x
2

x
0

x
1

x
2

a

б

y =

y
(
x
)

y = L
2
(
x
)

y = L
2
(
x
)

y = L
3
(
x
)

+

_

Рис. 6.5

Остаточный член общей формулы Сим
п
сона 6.16 имеет порядок
O
(
h
4
).




0. Формулы Ньютона

Котеса высших порядков

Полагая в формул
е

(6.6)
п

 , можно вычи
с
лить з
на
чения коэффициентов
A
i
:

A
0

=
,
A
1

=
A
2

=
,
A
3

=
.

и получить квадратурную формулу Ньютона




. (6.1
7
)

Формул
а 6.17 имеет погрешность того же порядка, что и формула Симпсона 6.15, т. е.
O
(
h
5
).

Приведем таблицу значений коэффицие
н
тов Котеса табл.
1
).

Таблица 1

п

Н
о

Н
1

H
2

H
з

H
4

H
5

H
6

Множи
тель для
п
о
лу
чения

А
i

1

1
/
2

1
/
2






h

2

1
/
6

4
/
6

1
/
6





2
h

3

1
/
8

3
/
8

3
/
8

1
/
8




3
h


Окончание табл.
1

п

Н
о

Н
1

H
2

H
з

H
4

H
5

H
6

Множи
тель для
п
о
лу
чения

А
i

4

7
/
90

32
/
90

12
/
90

32
/
90

7
/
90



4
h

5

19
/
288

75
/
288

50
/
288

50
/
288

75
/
288

19
/
288


5
h

6

41
/
840

216
/
840

27
/
840

272
/
840

27
/
840

216
/
840

41
/
840

6
h


Например, квадратурна
я формула Ньютона

Кот
е
са для
п

 5 имеет вид




.
(6.1
8
)

Квадратурные формулы с нечетным числом у
з
лов 
п

=
2, 4, 6 являются более удобными, т.
е. выгоднее при
менять формулу Симпсона 6.15, че
м формулу Ньютона 6.17, так как при
одном и том же порядке погрешности формула Нь
ю
тона требует больше узлов и больше
вычис
лений, чем формула Симпсона. Аналогично, формула для
п

=

4 лучше, чем форм
у
ла
для
п

 5, и т. д.




31.
Правило Рунге

оценки

погрешности квадратурных формул

Приведем сводку квадратурных формул с о
с
таточны
ми членами.

Формула трап
е
ций:



.




.


Рассмотрим на примере общей формулы Симпсон
а 6.16 правило Рунге для оценки
погрешности. Пусть подынтегральная функция четырежды непр
е
рывно диф
ференцируема.
Тогда формула остаточного члена им
е
ет вид

R
(
h
) =

(6.19)


где
ξ



некоторое число из о
т
резка [
а
,
b
].

Предположим, что производная
y
IV
(
x
 изменяется на этом отрезке медленно, и приближенно
можно записать остаточный член в виде
R
(
h
) =
М
h
4
, где
М



пост
о
янная. Пусть
S
h
, и
S
2
h
,
соответственно значения интеграла, полу
ченные по о
б
щей формуле Симпсона с ш
агом
h

и 2
h
.
Тог
да справедливы соотн
о
шения


=
,

=
.

Отсюда получим формулу для оценки погрешн
о
сти


(6.20)

Уточненное значен
ие интеграла по формуле Симпсона

вычисляется с учетом поправки




=

(6.21)

Для формулы трапеций также можно применить пра
вило Рунге. Так как фо
р
мула
остаточного члена общей формул
ы трапеций м
о
жет быть представлена в виде
R
(
h
) =
Mh
2
,
справедливы соотношения


=
,


=
.

Отсюда получим формулу для оценки п
о
грешности:

Mh
2

=

R
(
h
) =

.

Теперь для интеграла можно записать по правилу Рунге формулу трапеций с поправкой


=


(6.22)

Замечание.

Здесь необходимо подчеркнуть, что описанное правило

Рунге применимо
только тогда, когда выполняются указанные выше условия для функции
у
(
х
)


существование производной соответству
ю
щего порядка и ее ограниченность, точнее
говоря, возможность приближенного представл
е
ния по
грешности в виде
R
(
h
) =
Mh
4

для
фор
мулы Симпсона 
R
(
h
) =
Mh
2
для форм
у
лы трапеций, где
М



постоянная. Погрешность
представления остаточн
о
го члена в указанном виде считается достаточно малой. Эти
условия для конкре
т
ной функции могут не выполняться, тогда правило Рунге не будет
гарантир
о
ват
ь приемлемый резул
ь
тат.





2. Задача Коши и её решение методом Эйлера ломаных и м
етод Эйлера с
уточнением

Рассмотрим задачу Коши для обыкновенного диффе
ренциального уравнения пе
р
вого порядка

y
'
(
х
) =
f
(
x
,
у
),




(7.1)

y
(
x
0
) =
y
0


(7
.2)

Требуется найти функцию
у

=
у
(
х
, которая удовлет
воряет уравнению 7.1 на и
н
тервале 
х
0
,
X
 и
начально
му условию 7.2 в то
ч
ке
х
0
.

Приведем без доказательства теорему существования и единственности решения задачи К
о
ши.

Теорема
7.1.

Пусть в области

R
{(
x
,
у
),
|
х

-

х
0
|



а
,
|
у

-

у
0
|



b
} функция
f
(
x
,
у
 непр
е
рывна. Тогда на
некото
ром отрезке
|
х

-

x
0
|



d

существует решение уравнения 7.1, удовлетворяющее усл
о
вию 7.2.

Если в области
R

функция


f
(
x
,
у
 удовл
е
творяет усло
вию Липшица

|
f
(
x
,
у
1
)
-

f
(
x
,

y
2
)
|



k

|
y
1

-

у
2
|
,
k


� 0,

k


=
const
,

то указанное решение единс
т
венно.

Произведем разбиение отрезка [
х
0
,
X
] на
п

ча
с
тей:

x
i

=
x
0

+

ih
,
.


(7.3)

Найдем приближенные значения р
е
шения


у
(
х
)

в точ
ках

x
i

,
i

=

1,

2,..
.
,
n
.

1

Метод Эйлера
(
лом
а
ных
)

Рассмотрим уравнение 7.1 в точках

x
i
,

i

=

0, 1, ... ,
n

-

1

и

заменим производную
y
'
(
x
i
)

разностной
форм
у
лой


y
'
(
x
i
)

=

(?)


(7.4)

Тогда получим рекуррентную формулу
м
етода Эйле
ра

для вычисления прибл
и
женных
знач
е
ний
у
(
x
i

+

1
):

y
i
+ 1
=
y
i

+

h f
(
х
i
,
у
i
),

i

= 0, 1,...,
п

-

1.
(7.5)

Здесь через
y
i

обозначены приближенные зн
а
чения
y
(
x
i
, т. е.

у
i

=
у

(
x
i
),

i


= 1, 2, ... ,
п
.

На рис. 7.1 дана геометрич
еская иллюстрация метода Эйлера 7.5. Уравнение кас
а
тельной к
графику решения
у
(
х
 в точке 
х
0
,
у
0
 имеет вид

y

=
y
0

+
f
(
x
0
,
y
0
)(
x



x
0
)

(7.6)

так как
у
'
(
х
0
) =
f
(
x
0
,
y
0
. Интегральная кривая

у
(
х
 на от
резке [
х
0
,
х
1
] заменяется отре
з
ком касательной
(7.6)
, со
единяющей точку 
х
0
,
у
0
 с точкой 
х
1
,

у
1
, где


у
1

=
у
0

+
f
(
x
0
,
у
0
)(
х

-

х
0
 рис. 7.1. Точка 
х
1
,
у
1
)
уже не лежит на интегральной кривой
у

=
у
(
х
, удо
в
летв
о
ряющей на
чальному условию 7.2.

При
i

=

1 формула 7.5 дает точку 
х
2
,
у
2
, которая оп
ре
деляется с помощью кас
а
тельной

у

=
y
1
+
f
(
x
1
,

у
1
)(
х

-

x
1
, проведенной в точке 
х
1
,
у
1
 к инт
е
гральной кривой
у
(
х
),





Рис. 7.1

удовлетворяющей
уравнению 7.1 и начальному у
с
ловию

y
(
x
1
) =
y
1
.

Таким образом, с каждым шагом
i

метод Эйлера 7.5 дает точки 
x
i
,

y
i
, которые,
вообще говоря, удаляются от интегральной кривой, соответс
т
вующей точному
решению задачи Коши 7.1, 7.2. Вместо интеграл
ь
ной

кривой ме
тод Эйлера дает
ломаную, изображе
н
ную на рис. 7.1, по
этому его часто н
а
зывают
методом ломаных
.

Формулу 7.5 можно получить и другим способом. Рассмотрим разл
о
жение
искомого решения
у
(
х
 в ряд Тейлора в окрестности точки
х
0
:

у
(
х
) =
у
(
х
0
) +
у

(
х
0
)(
х

-

х
0
) +

+ ... .

Ограничившись двумя слагаемыми и учитывая, что
y
'
(
х
0
) =
f
(
x
0
,
у
0
, при
х

=
х
1

получим 7.5.

Т
акже
здесь
получ
ен

следующий результат


погреш
ность вычисления значения
у
1

есть величина порядка
O
(
h
2
, а погрешност
ь значения
у
п



вел
и
чина порядка
O
(
h
. Из
-
за
большой погрешности метод Эйлера прим
е
ня
ется редко.

2

Метод Эйлера с уточнением

Для повышения точности метода Эйлера применяют следующий прием. Сн
а
чала находят
приближенное зна
чение решения по методу Э
й
лера:

,

а затем уточняют его по формуле

.

Этот метод называется
методом Эйлера
-
Коши
, и ре
куррентные соотнош
е
ния

для его
реализации могут быть з
а
писаны в виде

, (7.7)

i

= 0, 1, ...,
n

-

1

Метод Эйлера
-
Коши имеет погрешность порядка
O
(
h
2
. Геометрическая и
л
люстрация
метода Эйлера
-
Коши пока
зана на рис. 7.2. Очередное прибл
и
жение метода Эйлера
-
Коши
соответствует точке пересечения диагоналей п
а
рал
лелограмма, построенного на двух звеньях
л
о
ман
ой метода Эйлера.



Метод Эйлера
-
Коши является одним из час
т
ных слу
чаев методов Рунге
-
Кутта.

. Методы Рунге
-
Кутта решения обыкновенных дифференциальных
уравнений

Рассмотрим наиболее популярный метод решения за
дачи Коши


метод Ру
н
ге
-
Кутта. Этот
метод

позволяет строить формулы расчета приближенного реш
е
ния прак
тически любого
п
о
рядка точности.

Выведем формулы метода Рунге
-
Кутта второго поряд
ка точности. Для этого решение
представим куском ряда Тейлора, отбр
а
сывая члены с порядком выше второго. Тогда
п
риближенное значение искомой функции в точке
x
i
, можно з
а
писать в виде:

y
1

=
у
(
х
0
) +
у

(
х
0
)(
х

-

х
0
) +

=
y
0

+
f
(
x
0
,
y
0
)
h

+
.
(7.8)

Вторую производную у
"
(
х
0
 можно выразить через производную функции
f
(
x
,
у
, однако в
методе Рунге
-
Кут
та вместо производной используют ра
з
ность


,

соответственно подбирая значения параметров

,
,
.
Тогда 7.8 можно п
е
реписать в
виде

y
1

=
y
0

+
h
[

f
(
x
0
,

y
0
) +
α
f
(
x
0

+
γ
h
,
y
0

+
δ
h
)],


(7.9)

где
α
,
β
,
γ

и
δ



некоторые пар
а
метры.

Рассматривая правую часть 7.9 как функцию аргу
мента
h
,
разложим ее по степ
е
ням
h
:

y
1

=
y
0

+
(
α

+
β
)
h
f
(
x
0
,

y
0
) +
α
h
2
[
γ

f
x
(
x
0
,
y
0
)

+

δ
f
(
x
0
,
y
0
)]
,

и выбере
м параметры
α
,
β
,
γ

и
δ

так, чтобы это разложе
ние было близко к 7.8. Отсюда
след
у
ет, что

α

+
β

= 1,
αγ


= 0,5,
αδ

=
0,5
f
(
х
0
,
y
0
).

С помощью этих уравнений выразим
β
,
γ

и
δ

через па
раметр
α
, получим

y
1

=
y
0

+
h
[
(1


α
)

f
(
x
0
,

y
0
) +
α
f
(
x
0

+
,
y
0

+
f
(
x
0
,
y
0
)
)
]
,

0
α



1. (7.10)

Теперь, если вместо
(
х
0
,

у
0
)

в 7.10 по
д
ставить 
х
1
,
у
1
),

получим формулу для вычисления
у
2



прибл
и
женного значения искомой функции в точке
х
2
.

В
общем случае метод Рунге
-
Кутта пр
и
меняется на произвольном разбиении отрезка [
х
0
,
X
]

на
п
частей, т. е. с переменным ш
а
гом

x
0
,
x
1
, ...
,
x
n
;
h
i

=
x
i

+ 1


x
i
,
x
n

=
X
.
(7.11)

Параметр
α

выбирают равным 1 или 0,5. Запишем
окончат
ельно
расчетные формулы
метода Рунге
-
Кутта второго порядка с пер
е
менным шагом для
α

= 1
:

y
i
+

1

=
y
i

+
h
i

f
(
x
i

+
,
y
i

+
f
(
x
i
,
y
i
)
)
,

(7.12)

i

= 0, 1, ...,
n



1
.

и

α  0,5:

y
i
+

1

=
y
i

+
[
f
(
x
i
,
y
i
)

+
f
(
x
i

+
h
i
,
y
i

+
h
i
f
(
x
i
,
y
i
))], (7.1
3
)

i

= 0, 1, ...,
n



1
.

Дадим геометрическую интерпретацию мет
о
дов Рун
ге
-
Кутта 7.12, 7.1.

На рис. 7.
а

иллюстрируется формула 7.12. Снача
ла по методу Эйлера в
ы
числ
я
ется
при
бли
женное значен
ие в точке
х
i

+
h
i
/2. В этой точке определяется касатель
ная к
интегральной кривой, параллельно которой через точку
(
x
i
,

y
i
)

проводится прямая до точки
пересечения с прямой

х
=
x
i

+

1
.

Ордината точки п
е
ресечения
y
i

+

1

при
нимается за
приб
лиженное значение искомого реше
ния в то
ч
ке


x
i

+

1
.

Рис. 7.,
б
интерпретирует формулу 7.1. Методом
Эй
лера вычисляется приближенное
значение
y
i

+
h
i
f
(
х
i
,
y
i
)

в
точке
х
i

+

1
.

В этой точке тангенс угла н
а
клона касатель
ной к
интегральной кривой равен выр
аж
е
нию


f
(
x
i

+
h
i
,

y
i

+
h
i
f
(
х
i
,
y
i
)).

Через точку 
x
i
,
у
i
 проводится
прямая, тан
генс угла наклона которой опр
е
деляется как полусумма угловых коэффициентов
касательных, проведенных ч
е
рез точки
(
x
i
,

y
i
 и 
x
i

+
1
,

y
i

+
h
i
f
(
x
i
,
у
i
)).

За приближенное зна
чен
ие
искомого р
е
шения в точке
х
i

+1

принимается орди
ната точки пересечения этой прямой с
пр
я
мой
х

=
x
i
+
1
.

Отметим, что метод
7.1 совпадает с мет
о
дом Эйле
ра
-
Коши 7.7.






Рис. 7.

Теоре
ма
7
.
2
.

Если правая часть
f
(
х
,
у
)

уравнения 7.1 непрерывна и огран
и
чена вместе со
своими производными до второго порядка включительно, то р
е
шение, получен
ное по
формулам 7.12 и 7.1, равномерно сходится к реш
е
нию задачи 7.1 и 7.2 с
п
о
грешностью
O
(
max

h
i
2
).

Приведем наиболее употребительные формулы метода Рунге
-
Кутта, форм
у
лы четвертого
порядка то
ч
ности:



(7.14)






4. Метод Рунге
-
Кутта с автоматическим выбором шага. П
равило Рунге
оценки погрешности

Для метода Рунге
-
Кутта применимо пр
а
вило Рунге для оценки погрешности. Пусть
у
(
х
;

h
)
-

приближенное значение решения в точке
х
, полученное по фо
р
мулам 7.12, 7.1 или 7.14
с шагом
h
, а

р

-

порядок точнос
ти соответству
ю
щей формулы. Тогда погрешность
R
(
h
)
значения
у
(
х
;
h
 можно оценить, испол
ь
зуя приближен
ное значение
у
(
х
; 2
h
 решения в точке
х
, полученное с ш
а
гом 2
h
:

R
(
h
) =

(7.15)

где
р

= 2

для
форм
ул 7.12 и 7.1 и
р

 4 для
(7.14).

Уточненное решение запишем в виде



(7.16)

В алгоритмах с автоматическим выбором шага предва
рительно задают п
о
грешность в виде
положительного параметра
ε
, и на каждом этапе
вычисл
е
ния следующего значения
y
i

+

1

подбирают такой шаг
h
, при котором вы
полняется неравенство


(7.17)


5. Задача Коши для системы дифференциальных уравнений и её
решение

Метод Рунге
-
Кутта применим и к

задаче Коши для системы
m

диффере
н
циальных
уравнений первого поряд
к
a

с
m

неизвестн
ы
ми функциями


x



(
x
0
,
X
)
(7.18)

y
1
(
x
0
) =
y
1,0
,
y
2
(
x
0
) =
y
2,0
,…
,
y
m
(
x
0
) =
y
m
,0

(7.19)

Приведем для задачи 7.18, 7
.19 расче
т
ные форму
лы метода Рунге
-
Кутта четвертого
порядка. Пусть требу
ется найти систему
т

функций
у
1
(
х
),
у
2
(
х
),...
,
у
m
(
х
, удовлетворяющих
в интервале 
х
0
,
X
 дифф
е
ренциаль
ным уравнениям 7.18, а в точке
х
0

-

начальным усло
-
виям 7.19. Пре
д
положим
, что отрезок [
х
0
,
X
] разбит на
N

частей:

x
i

=
х
0

+
ih
,

h

=
.


Тогда каждую
l
-
ю функцию
у
l
(
х
 можно приближенно вычислять в то
ч
ках
х
i

+1

по
формулам Ру
н
ге
-
Кутта

K
l
,1

=
f
l
(
x
i
,

y
1,
i
,
y
2
,
i
,…

,
y
m
,
i
)
,

i

= 1, 2, ... ,
m
,

K
l
,2

=
f
l
(
x
i

+
,

y
1,
i

+
,
y
2
,
i
+
,…

,
y
m
,
i

+
)
,



i

= 1, 2, ... ,
m
;

K
l
,3

=
f
l
(
x
i

+
,

y
1,
i

+
,
y
2
,
i
+
,…

,
y
m
,
i

+
)
,



i

= 1, 2, ... ,
m
;

K
l
,4

=
f
l
(
x
i

+
h
,

y
1,
i

+
,
y
2
,
i
+
,…

,
y
m
,
i

+
)
,




i

= 1, 2, ... ,
m
;

y
l
,
i
+1

=
y
l
,i

+
(
+

+

+
)
,


(7.20)


i

= 1, 2, ... ,
m
;


Здесь через
у
l
,
i

обозначено

приближенное значение функции
у
l
(
х
 в точке
х
i
.

Необходимо обратить

внимание на порядок вычислений по фор
мулам 7.20. На
каждом шаге сначала вычисляются ко
эффициенты
K
l
,
i

в следу
ю
щем порядке:



и лишь затем приближе
нные значения фун
к
ций
y
1,
i
+1
,
y
2,
i
+1
,…,
y
m
,
i
+1
.


6. Задачи Коши для дифференциальных уравнений высших порядков

Задача Коши для дифференциального ура
в
нения
n
-
го порядка

y
(
n
)

=

f
(
x
,
у
,
у
', ... ,
y
(
n
-
1)
),

x



(
х
0
,
X
),

(7.21)

y
(
x
0
) =
y
0
,
y
'
(
x
0
) =
y
1,0
,…,
y
(
п
-
1)
(
x
0
) =
y
п

-

1,0
.


(7.22)

легко сводится к задаче Коши для системы дифференци
альных уравнений первого
порядка с помощью замены пер
е
менных

z
0

=

y
,

z
1

=
y
'
,...,

z
n

-

1

=
y
(
n



1
)
.


(7.23)

Учитывая 7.2, из ур
авнения 7.21 п
о
лучим систе
му дифференциальных уравн
е
ний


(7.24)

Начальные условия 7.22 для функций
z
l
, п
е
реписыва
ются в виде

z
0
(
x
0
) =
y
0
,
z
1
(
x
0
) =
y
1,0
,…,
z
n
-
1
(
x
0
) =
y
п

-

1,0
.


(7.25)

З
апишем для полученной системы метод Рунге
-
Кутта:

z
l
,
i
+1

=
z
l
,
i

+
(
+

+

+
), (7.26)


i

= 0, 1, 2, ...
,
N
;
l

=
0, 1,…
,
n



1.

Для вычисления коэффициентов
,
,
,

имеем следующие формулы:

K
0,1

=
z
1,
i
,

K
1
,1

=
z
2
,
i
,

…………

K
n
-
1
,1

=
f
(
x
i
,
z
0
,
i
,
z
1,
i
,…,
z
n
-
1,
i
),

K
0,2

=
z
1,
i

+
,

K
1
,2

=
z
2
,
i

+
,

…………………

K
n
-
1
,2

=
f
(
x
i

+
,
z
0
,
i

+
,
z
1,
i

+

,…,
z
n
-
1,
i

+
),

K
0,
3

=
z
1,
i

+
,

K
1
,
3

=
z
2
,
i

+
,

…………………

K
n
-
1
,3

=
f
(
x
i

+
,
z
0
,
i

+
,
z
1,
i

+

,…,
z
n
-
1,
i

+
),

K
0,
4

=
z
1,
i

+
,

K
1
,
4

=
z
2
,
i

+
,

K
n
-
1
,4

=
f
(
x
i

+
,
z
0
,
i

+
,
z
1,
i

+

,…,
z
n
-
1,
i

+
).


7. Краевые задачи для обыкновенных дифференциальных уравнений.
Метод прогонки

Пр
имерами краевых задач для обыкновенных диффе
ренциальных уравн
е
ний являются
сл
е
дующие:

1 Найти функцию
u
(
х
, которая удовлетворяет на от
резке [
а
,
b
]

уравн
е
нию

u
"
(
х
) =
-
f
(
х
)
,

(7.27)

а

на концах отрезка ус
ловиям

u
(
а
)

=

u

(
b
) =
0
.


(7.28)


Задача 7.27, 7.28 имеет следующее физическое со
держание. Между то
ч
ками
х

=
а

и
х

=
b

натянута упру
гая струна, находящаяся под действием внешней изгиба
ющей нагрузки.
f
(
x
)
-

величина нагрузки, а
u
(
х
)
-

прогиб стр
у
ны в безразмерных единицах.

2)

Двухточечная краевая задача для линейного диффе
ренциального ура
в
нения второго
п
о
рядка:

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) =
f
(
x
),


(7.29)


α
0
u
(
a
) +
α
1
u
'
(
a
) =
A
,

β
0
u
(
b
) +
β
1
u
'
(
b
) =
B
.


(7.30)

Метод прогонки

Рассмотрим частный случай задачи 7.29, 7.0:

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) =
f
(
x
), (7.31)

u
(0)
=
а
,
u
(
Х
) =
b
.

(7.32)

Введем сетку
x
i

=
ih
,

i

= 1, 2, ... ,
N
;

h

=
X
/
N
. Обозна
чим через
u
i

прибл
и
женные значения
u
(
x
 в узл
ах сетки. Рассмотрим уравнение 7.1 во вну
т
ренних узлах
х
i

i

=

1,

2,
...,

N

-


заменим производную второго поряд
ка ра
з
ностной формулой 5.11


(7.33)

Тогда из 7.1, 7.2 получим для определения
u
i

си
стему линейных ура
в
нений



(7.34)

Система 7.4 при
р
(
x
)


0 имеет решение. Система 7.4 представляет с
о
бой систему
линейных уравнений с трехдиагональной матрицей


(7.35)

и для ее решения применим мет
од прогонки, который фактически является методом
исключения неизвестных Гау
с
са.

1
.
Прямой ход прогонки
. Запишем первое уравнение 7.5 в виде

.
(7.36)

Подставив 7.6 во второе уравнение 7.5 и упрос
тив выражения, ув
и
дим, что мо
жно
для
u
i

пол
у
чить фор
мулы


(7.37)

Из уравнения 7.5, учитывая 7.7 при
i

=
N

-

2, пол
у
чим


(7.38)

2.
Обратный ход прогонки.

После вычис
ления прогон
очных коэффицие
н
тов можно
на
й
ти

значения решения
за
дачи. Формула 7.8 дает значение

u
N

-

1
:

u
N

-

1

=
β
N

-

1
(7.39)



Остальные значения вычисляем в обратном поря
д
ке
:




(7.40)


3
8. Краевые задачи для обыкновенных дифференциальных
уравнений.
Метод стрельбы пристрелки

Рассмотрим метод пристрелки на примере решения линейной краевой з
а
дачи 7.1,
(7.32):

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) =
f
(
x
),


u
(0)
=
а
,
u
(
Х
) =
b
.



Если известны частное решение
u
(
х
 нео
д
нородного уравнения

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) =
f
(
x
),

и два линейно независи
мых реш
е
ния
и
1
(
х
),
u
2
(
х
 однородного уравнения

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) = 0,

то общее решение неоднородного диффе
ренциального

ура
в
нения

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) =
f
(
x
),


можно за
писать в виде

u
(
x
) +
С
1
u
1
(
х
) +
C
2
u
2
(
x
).

Постоянные
С
1
,
С
2

можно определить из краевых ус
ловий 7.2.

В методе пристрелки используется сл
е
дующий способ.


Сначала находят частное решение
u
(
х
 н
еоднородного уравнения

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) =
f
(
x
),

удовлетворяющее усло
вию
u
(0) =
а
, и частное решение
u
1
(
х
 однородного урав
нения

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) = 0,

удовлетв
о
ряющее условию
u
(0) = 0.


Затем общее решение
u
(
х
 неоднородного
уравнения

u
"
(
x
) +
p
(
x
)
u
'
(
x
) +
q
(
x
)
u
(
x
) =
f
(
x
),

удовлетворяющее условию
u
(0) =
а
, записывают в виде
u
(
х
) +
С
u
1
(
х
. Остае
т
ся найти
постоян
ную
С
1

из условия
u
(
x
) +
C
u
1
(
x
) =
b
.

Приведем сеточный аналог метода пристрелки. Пусть краевая задача пр
и
ведена к
сист
еме уравнений 7.4


(7.41)

Найдем частные решения неоднородной 7.42 и одно
родной 7.4 си
с
тем уравнений:


(7.42)


(7.43)

Выбирая прои
звольные значения для

при этом

должно быть

, находим из
7.42 и 7.4 формулы для вычисления частных реш
е
ний:


(7.44)




(7.45)


Найдем
С

из условия
и зап
и
шем решение:


(7.46)


(7.47)

Алгоритм метода пристрелки заключается в том, что
б
ы выбрать шаг
h

и выполнить
последовател
ь
но

вычисле
ния по формулам 7.44

-

(7.47).


9. Дифференциальные у
равнения в частных производных

в
олновое уравнение уравнение колебаний, уравнение
теплопроводности и уравнение
Пуассона

Уравнение


F
(
х
1
,

х
2
,...,

х
m
,

u
,

,
...,
,
,

,...) = 0,

с
в
язывающее неизвестную функцию
u
(
х
1
,

х
2
,...,

х
m
),
не
зависимые переменные
х
1
,

х
2
,...,

х
m


и частные произ
в
одные
неизвестной функции, называется
ура
в
нением в частных
производных.
Порядок
п

старшей производной называется порядком дифференциального
ура
в
нения.

Уравнение в частных производных называется линей
ным, если оно л
и
нейно
относ
и
тельно неизвестной функ
ции и е
ё

частных производных.

Приведем основные уравнения математ
и
ческой физ
и
ки, которые являются
линейными уравнениями в част
ных производных второго порядка.

Волновое уравнение


(
уравнение колеб
а
ний
)


(8.1)

описывает различные виды волн



звуковые, упругие,
эл
ектромагнитные и другие
к
олебательные процессы. Функция
u

=
u
(
х
,

у
,

z
,

t
 зависит от пр
о
странственных
переменных
х
,

у
,
z

и времени
t
.



2
.

Процесс распространения тепла в одноро
д
ном изот
ропном теле

описыв
а
ется
уравнением теплопроводности





(8.2)


Уравнение теплопроводности

в общем виде записыва
ется так:

c
(
х
,

у
,

z
,

t
)

=
div
[
k
(
u
,

х
,

у
,

z
,

t
)

grad

u
] +
q
(
х
,

у
,

z
,

t
)
,

grad

u

=
,

div(
X
(
x
,

y
,

z
),
Y
(
x
,
y
,
z
),

Z
(
x
,

y
,

z
)) =

Здесь
u

=
u
(
х
,

у
,

z
,

t
)


температура в точке 
х
,

у
,

z
 в момент времени
t
,

c
(
u
,

х
,

у
,

z
,

t
)



теплоемкость в точке 
х
,

у
,

z
 в момент времени
t
,
k
(
u
,

х
,

у
,

z
,

t
)


коэффици
-
ент теплопроводности в то
ч
ке 
х
,

у
,

z
 в момент времени
t
,

q
(
u
,

х
,

у
,

z
,

t
)


плотность источников

тепла.


3
.

Установившееся тепловое состояние в однородном изотропном теле

оп
и
сывается
уравнением Пуассона

=
-
f
(
x
,
y
,
z
). (8.3)

Установившееся тепловое состояние в однородном изотропном теле при
отсутствии
источн
и
ков тепла внут
ри тела описывается
уравнением Лапласа

= 0.


(8.4)

В уравнениях 8. и 8.4 функция
u

=
и
(
х
,

у
,

z
 зави
сит только от пр
о
странственных
п
е
ременных
х
,

у
,

z
.

Согласно
классификации уравнений второго порядка уравнение 8.1 отн
о
сится к
гиперболическому типу
, (8.2)


к
параболическому типу
, а 8., 8.4


к
эллип
-
тическому типу
.

В конкретной постановке задачи математической фи
зики необходимо на
й
ти решение
одного из ура
внений 8.1



8.4, удовлетворяющее дополнител
ь
ным 
началь
ным

и
гр
а
ничным
 условиям.

Начальные условия

задаются с уравнениями 8.1, 8.2 и обычно имеют вид

u
(
x
,

y
,

z
,

t
0
)

=
u
0
(
x
,

y
,

z
),

(8.5)



=
u
1
(
x
,

y
,

z
).
(8.6)

При этом для уравнения 8.1 искомое решение
u

=
u
(
х
,

у
,

z
,

t
 должно удовлетворять в
начальный момент времени
t
0

условиям 8.5 и 8.6, а для уравнения 8.2


одному из
у
с
ловий 8.5, 8.6.

Граничные условия

для уравнен
ий 8.1, 8.2: искомое решение

u

=
u
(
х
,

у
,

z
,

t
 должно удовлетворять на грани
це тела или среды одному из условий

, (8.7)


,

(8.8)

где



производна
я

по направлению норм
а
ли к грани
це
G

тела.

Для уравнений Пуассона или Лапласа з
а
даются толь
ко
граничные

условия



(8.9)

или



(8.10)

Задача для уравнения Лапласа 8.5 с гр
а
ничным ус
ло
вием 8.9 называется
задачей
Дирихле
, а с условием 8.10


задачей Не
й
мана
.

Область
Ω
, описывающая тело, может быть задана в трехмерном, двуме
р
ном или
одномерном пространствах.
В

первом случае функция
u

=
u
(
х
,

у
,

z
,

t
 зависит от трех
пр
о
странственных пе
ременных 
х
,

у
,

z
 и времени
t
, во втором


от двух
пространственных переменных 
х
,
у
 и времени
t
:
u

=
u
(
х
,

у
,

t
, а в третьем случае


от
перемен
ных
х

и
t
:
и  и
(
х
,

t
)
.

Р
ассмотрим
разностный метод

решения задач для уравнений в частных
производных, котор
ый описывает
ся ниже на примерах основных задач мат
е
матической
физики.


40.
Разностный метод для уравнения колебаний
.
Явная схема

Рассмотрим задачу о малых колебаниях натянутой стру
ны с распределе
н
ной по длине
нагрузкой
f
(
x
,
t
 рис. 8.1:

u
tt

=
c
2
u
xx

+
f
(
x
,

t
),

0



x



a
,

0



t



T
, (8.11)

u
(
x
, 0) =
μ
(
x
),
u
t
(
x
, 0) =
μ
0
(x), 0
x


a
, (8.12)

u
(0,
t
) =
μ
1
(
t
),
u
2
(
a
,
t
) =
μ
2
(
t
), 0



t



T
.
(8.13)



Рис. 8.1

Струна совершает
плоские

колебания
, т. е. точки стру
ны перемещаются параллель
но
плоскости
t

= 0.

Функция
u
(
х
,
t
 выражает смещение точки
х

струны в момент времени
t

от
прямолине
й
ной формы.

Начальные условия

8.12 означают следующее. Фор
ма струны в начал
ь
ный момент
времени
t

 0 выражает
ся функцией
μ
(
x
)
. Скорость перемещения точк
и
х

стру
ны в
момент времени
t

 0 равна значению функции
μ
0
(
х
).

Краевые условия

8.1 говорят о том, что левый конец струны с течением времени
с
о
вершает смещение
μ
1
(
t
),
a

правый конец


смещение
μ
2
(
t
).

Если концы струны закреплены, то
μ
1
(
t
) =
μ
2
(
t
) = 0.

При этом

предполагае
тся
, что начальные условия 8.12 и кра
евые условия 8.1
должны быть согласованы между со
бой в угловых точках, т. е. выпо
л
нены условия
u
(0,
0) =
μ
(0)

=
μ
1
(
0
),
u
(
а
, 0)
=

μ
(
a
)

=
μ
2
(
a
).


На рис. 8.1 представлен случай, когда

u
(0, 0) =

μ
(0)

=
μ
1
(
0
)

= 0,

u
(
а
, 0)
=

μ
(
a
)

=
μ
2
(
a
) = 0.

Введем сеточную область рис. 8.2, а. В прямоуголь
ной области
0


x



a
,

0


t



T


з
а
дадим точки

(
x
i
,
t
k
),
x
i

=
ih
,
i

=
l
,...,
N
;

t
k

=
k
τ
,

k

= 0, 1, ... ,
M
;
τ

=



(8.14)

Рассмотрим уравнение 8.11 в точках
(
x
i
,
t
k
)
,
i

= 1,...
,
N

-

1,


k

= 1,...
,
М

-

1,

и заменим
производные разно
стными формулами



u
tt
(
x
i
,

t
k
)

=




(8.15)

u
xx
(
x
i
,
t
k
)

=


Обозначим через
u
i
,
k

приближенные знач
е
ния искомой функции в точках 
x
i
,
t
k
. Тогда
из уравнения 8.11 п
олу
чим
разностное уравнение

(
разнос
т
ную схему
, которое
аппро
к
симирует

уравнение 8.11 с порядком
О
(
h
2

+

τ
2
):



(8.16)

i

= 1,...,
N

-

1;
k

= 1,...,
M

-

1.


На рис
. 8.2,
б

изображен шаблон крестª разностного уравнения 8.16. Разностное
уравнение 8.16 связывает значения неизвестной функции на трех слоях 
k



1,
k
,
k

+ 1)





На слое
k

 0 заданы начальные условия 8.12, из ко
торых следует, что

u
i
,

0

=
μ
(
x
i
)
,
i

= 1, 2, ... ,

N

-

1.

(8.17)

Чтобы найти значения неизвестной функции на слое
k

 1, используем у
с
ловие для
производной
u
t
(
x
,

0 из 8.12. Для этого п
о
строим разложение в ряд Тейлора:

u
(
x
i
,

t
1
) =
u
(
x
i
,

τ
) =
u
(
x
i
,

0) +
u
t
(
х
i
,

0)
τ

+

u
tt
(
x
i
,

0)
τ
2
+
O
(
τ
3
).


(8.18)

Из уравнения 8.11, учитывая первое условие в 8.12, выразим вторую производную

u
tt
(
x
i
,

0)

=
c
2
u
xx
(
x
i
,

0) +

f
(
x
i
,

0)

=
с
2
μ''
(
х
i
) +
f
(
х
i
,

0).



(8.19)

Теперь, учитывая условие
u
t
(
x
,

0) =
μ
0
(
x
)
в 8.12, из

8.18, 8.19 выводим формулу
для вычисления значений

функции на первом слое:

u
(
x
i
,
t
1
) =
μ
(
x
i
) +
μ
0
(
x
i
)
τ

+
(
с
2
μ''
(
х
i
)

+
f
(
х
i
,

0)) +
О
(
τ
3
). (8.20)

С учетом 8.1 окончательно получим для приближен
ных з
начений и
с
комой функции
на первом слое формулы

u
i
,

1

=
μ
(
x
i
) +
μ
0
(
x
i
)
τ

+

i

= 1, 2,...,
N

-

1;
u
0,

1

=
μ
1
(
t
1
)
,
u
N
,

1

=
μ
2
(
t
1
).

(8.21)

Учитывая граничные условия 8.1,

из 8.16 выводим формулы для в
ы
числения
знач
ений на слоях
k

= 2,...
,
М
:

u
0,

k

=
μ
1
(
t
k
),
u
N
,

k

=
μ
2
(
t
k
)
.

u
i
,

k

=
2
u
i
,

k

-

1

-

u
i
,

k

-

2

+

i

= 1, 2,...,
N

-

1.
(8.22)

Таким образом
,

получ
ены

явные формулы 8.17, 8.21, 8.22 ре
шения разностной
з
а
дач
и.

Разностная схема называется
устойчивой
, если она имеет единственное решение и
м
а
лым изменениям исход
ных данных отвечают малые изменения решения.

Известен следующий факт
: для выполнения условия
устойчивости

ра
з
ностной
сх
е
мы 8.16 необходимо и достаточн
о, чтобы выполнялось
усло
вие Кура
н
та


с
τ


h
.

Г
овор
ят
, что решение разностной схемы 8.16
сходится

к решению и
с
ходной задачи
(8
.
11)



8.1, если
в
ыполняется условие

при
h



0
и
τ



0
.

Если разностная схема устойчива и аппроксим
ирует исходную задачу, решение
разностной схемы
сходится

к решению з
а
дачи.

Сформулируем алгоритм решения задачи о колебани
ях струны 8.11


8.1 с
пом
о
щью явной разностной схе
мы 8.16.

1.

Построить сеточную область

(
x
i
,
t
k
),
x
i

=
ih
,
i

= 0,
l
,...,
N
;

t
k

=
k
τ
,

k

= 0, 1, ... ,
M
;
τ

=


в
ыбирая шаги
h

и

τ

так, чт
о
бы выполнялось условие ус
тойчивости 
условие Куранта
)
сτ


h
.

2.

Вычислить значения
u
i
,
k

искомой фун
к
ции для
k

=

0,1:


u
i
,

0

=
μ
(
x
i
)
,

i

=
0, 1, ... ,
N
,


u
i
,

1

=
μ
(
x
i
) +
μ
0
(
x
i
)
τ

+

i

= 1,...,
N

-

1;
u
0,

1

=
μ
1
(
t
1
)
,
u
N
,

1

=
μ
2
(
t
1
).

(8.23)

3.
Вычислить

значения

u
i
,
k


для

k

= 2,...,
M
:


u
0,

k

=
μ
1
(
t
k
),
u
N
,

k

=
μ
2
(
t
k
)
.


u
i
,

k

= 2
u
i
,

k

-

1

-

u
i
,

k

-

2

+

i

= 1, 2,...,
N

-

1.

(8.24)

41. Уравнение колебаний струны. Неявная схема

Построим неявную схему для уравнения колебаний
с
тру
н
ы 8.11:



=



(8.25)

i

= 1,...,
N



1;
k

= 1,...,
M
;
.

Для устойчивости схемы 8.25 параметры
с
,

h
,

τ
,

σ

д
олжны удовлетв
о
рять условию:


Если 1/4


σ



1/2, схема 8.25 безусловно устойч
ива.

Шаблон схемы 8.25 изображен на рис. 8.4



Значения решения на нулевом и первом слоях вычис
ляют по формулам 8.2. На
каждом
k
-
м слое 
k

= 2, 3, ... ,
М
 решают методом прогонки сист
е
му уравнений
относительно
u
i
,
k
+1

i

= 1, 2, ... ,
N

-

1:



(8.26)



Алгоритм решения неявной разностной схемы для

уравнения кол
е
баний
струны
:

0. Построить сеточную область, выбирая шаги
h
,
τ

и па
раметр
σ

так, чтобы
выполн
я
лось
условие устойчивости


Если 1/4


σ



1/2, можно выбирать
h
,
τ

произвольно.

1.

Вычислить значения
u
i
,
k

искомой функции для
k

= 0, 1:


u
i
,0

=
μ
(
x
i
),
i

= 0, 1,...,
N
; (8.27)


u
i
,1

=
μ
(
x
i
) +

μ
0
(
x
i
)
τ

+
, (8.28)

i

= 1,...,
N



1;
u
0,

1

=
μ
1
(
t
1
),
u
N
,

1

=
μ
2
(
t
2
)

2.
Значения
u
i
,
k
+1

для каждого
k

 0 находим методом прогонки.

2.1. Вычислим правые части 8.26:





(8.29)

i

= 1,...,
N



1
.

2.2.
Вычислим прогоночные коэффициенты;



(8.30)



(8.31)


. (8.32)

2.. Вычислим решение
u
i
,
k
+1
:


u
0,
k

+ 1

=
μ
1
(
t
k

+ 1
),
u
N
,
k

+ 1

=
μ
2
(
t
k

+ 1
),


u
N



1
,
k

+ 1

=
β
N



1
, (8.33)


u
i
,
k

+ 1

=
α
i
u
i
+ 1
,
k

+ 1

+
β
i
,


i

=
N



2,
N



3,..., 1. (8.34)



42. Разностный метод для уравнения колебаний мембраны

Рассмотрим задачу для уравнения к
олебаний однород
ной прямоугольной мембраны:

u
tt

=
c
2
(
u
xx

+
u
yy
) +
f
(
x
,
у
,

t
),


(
8.35)

0
x


a
, 0
y


b
, 0
t


T
.

u
(
x
,
у
, 0) =
μ
(
x
,
у
), 0
х


a
, 0
у


b
. (8.36)

u
t
(
x
,
у
, 0) =
μ
0
(
x
,
у
), 0
х


a
, 0
у


b
.

(8.37)

u
(0,
у
,

t
) =
μ
1
(
y
,
t
),
u
(
a
,
у
,

t
) =
μ
2
(
y
,
t
)

0


y


b
, 0


t



T
.

u
(
x
,
0
,

t
) =
μ
3
(
x
,
t
),
u
(
x
,
b
,

t
) =
μ
4
(
x
,
t
)

(8.38)

0


x


a
, 0


t



T
.

Введем сеточную о
б
ласть:

(
x
i
,
y
j
,
t
k
),
x
i

=
ih
x
,
i

= 0, ...,
N
x
,
h
x

=
,

y
j

=
jh
y
,
j

= 0, ...,
N
y
,
h
y

=
,

t
k

=

,
k

= 0, ...,

M
,
τ
=
.

Обозначим

u
i,j,k

=
u
(
x
i
,
у
j
,
t
k
).
Заменяя производные Разностными формул
а
ми для
уравнения 8.5, получим
разностное уравнение с порядком аппро
к
симации

O
(
):



+

+

(8.39)

i

= 1, ...,
N
x



1,
j

= 1, ...,
N
y



1,
k

= 1, ...,
M



1.

Таким образом, получена явная

разностная схема, аналогичная явной сх
е
ме для
уравнения колебаний струны.


4. Вычисление определенного интеграла
методом Монте
-
Карло

В отличие от традиционных численных методов реше
ния задач, закл
ю
чающихся в
разработке алгоритма пост
роения решения,

для исследования некоторых классов задач
оказывается более целесообразным модел
и
рование их сущности с использованием
законов больших чисел теории вероятн
о
стей
. Оценки
f
1
,
f
2
, ... ,
f
n

искомой ве
личины
f

в этом

случае строят на осн
о
вании статистичес
кой
обработки материала, полученного в
результате мн
о
гократных случайных испытаний. Основным требо
ванием при этом
явл
я
ется сходимость по вероятности рассматриваемой случайной величины к иском
о
му
р
е
ше
нию задачи:


(9.1)

Здесь
Р



вероятн
ость,
ε



с
коль угодно малое поло
жительное число.

В отличие от численных методов, описанных ранее, в данном случае в
ы
числительный
процесс является
неде
терминированным
. Такой подход к решению вычисли
тельных
задач получил общее название
метода

Мо
н
те
-
Карло
. При практической реализации
данного подхода случайные величины заменяют сериями, так называемых
псевдослучайных величин
, генериру
е
мых соответствую
щими стандар
т
ными
программами
.

Вычисление определенного интегр
а
ла


Пусть задана непрерывная случайная ве
личина
ξ

(
кси
, с плотностью вер
о
ятности
р
(
х
),
значения которой распре
делены на и
н
тервале 
а
,
b
)

Плотность вероятности
р
(
х
 обладает сл
е
дующими свойствами:



1)
p
(
x
)


0;

2)
.

(
9.2)


Тогда математическое ожидание случа
й
ной величины
ξ

равно интегралу

M
(
ξ
) =
.

(
χ

-

хи

Для функции
f
(
ξ
, аргументом которой я
в
ляется слу
чайная величина
ξ
, т. е. для
случайной функции, матема
тическое ожидание ра
в
но

M
(
f
(
(
ξ
)
) =
.

(9.3)

Отсюда следует, что любой интеграл вида

где функция
р
(
х
 обладает свойствами 9.2, можно при
нять за математич
е
ское ожидание
некоторой сл
у
чайной функции
f
(
ξ
).

Но математическое ожид
ание случайной величины
f
(
ξ
 можно прибл
и
женно
определить с помощью статистичес
ких оценок, т. е. на основе выбо
р
ки {
ξ
1
,
ξ
2
,...,
ξ
N
}

объема
N

как среднее арифметическое зн
а
чений {
f
(
ξ
1
)
,
f
(
ξ
2
)
,...,
f
(
ξ
N
)
}
п
оэтому интеграл
(9.
3
 можно вычислить прибли
женно
по фо
р
муле





(9.4)

Теоретическим основанием для такого перехода явля
ется
центральная предельная
теорема

теории вероятно
стей, которая в упрощенной формул
и
ровке утверждает
сл
е
дующее.

С
реднее арифметическое
N

испытаний случайной ве
личины
ξ


ξ
N

=

является случайной величиной с тем же м
а
тематическим ожиданием


M
(
ξ
N
)
=

М
(
ξ
)


и
дисперсией,

р
авной

D
(
ξ
N
) =

и

при
N





закон распределения слу
чайной величины
ξ
N

стремится к но
р
мальному
зак
о
ну.

Очевидно, что чем больше
N
, тем меньше дисперсия
D
(
ξ
N
) =
. В
е
личину
погрешности формулы 9.4

можно
оценить по вероятности
. Н
а
пример, при боль
ших
N

можно утверждать, что с вероятн
о
стью 0,997 ошибка не превосходит величину 
σ

=
3


прави
ло трех сигмª для нормально распределенной случайной величины.

Можно считать, что погрешность формулы 9.4 есть

величи
на порядка

O
, но
для п
о
вышения точности
в

данном случае
нельзя применять правило Рунге
-
Ромберга
.

Приведем другой способ статистической оценки для одномерного инт
е
грала


.

Для этого вспомним его геометрический
смысл. Предпо
ложим, что фун
к
ция
f
(
x
)
непрерывна и
положительна на отрезке [
a
,
b
].

Тогда интеграл равен площади
криволинейной трапеции, ограниченной графиком функции
f
(
x
, осью аб
с
цисс и
прямыми
х

=
а

и
х

=
b

рис. 9.1.

Рассмотрим две случайные величины
ξ



равномерно распределенную на отрезке [
а
,
b
]

и
η

эта




равномерно рас
пределенную на о
т
резке [0,
f
max
]

где


f
max

=
.



Рис. 9.1

Вероятность попадания случайной точки 
ξ
,
η
 в криво
линейную т
рап
е
цию равна
отношению площади трапеции к площади пр
я
моугольника

{(
х
,
у
),
a



х



b
, 0


у



f
max

}:


(9.5)

Проведем серию из
N

испытаний и пол
у
чим
N

случай
ных точек 
х
,
у
, принадлежащих
прямоугольной облас
ти {
a



х



b
, 0


у



f
max

}. Подсчит
а
ем число
N
f

точек, удов
-
летворяющих условию
у



f
(
x
. Тогда вероятность поп
а
да
ния случайной точки 
ξ
,
η
 в
криволинейную трапецию приближенно ра
в
на относительной частоте попадания в

крив
о
линейную трапецию, т. е.
р




и интеграл при
ближенно вычисляе
т
ся по
формуле




(9.6)

Другой, более простой способ вычисления интеграла заключается в сл
е
дующем.

Проведем серию из
N

испытаний случай
ной величи
ны
,
равномерно ра
с
пределенной
на отрезке [
а
,
b
]
, и полу
чим
N

случайных чисел
x
i
, принадлеж
а
щих отрезку [
а
,
b
].
Вычислим инт
е
грал по формуле




(9.7)

Отметим, что в 9.7 подынтеграл
ьная функция может принимать полож
и
тельные и
отрицательные зн
а
чения, тогда как формула 9.6 применима только для неотрица
-
тельной функции
f
(
x
).

В общем случае, когда пределы интегрирования могут быть бесконечн
ы
ми,
необх
о
димо преобразовать интеграл к виду



(9.8)

и применить формулу 9.7.


44. Вычисление кратных интегралов
методом Монте
-
Карло

Метод Монте
-
Карло для вычисления о
д
номерных ин
тегралов обычно не
применяется, так как для получения высокой точности более удо
бны ква
д
ратурные
формулы. Этот метод оказывается более эффективным при в
ы
чис
лении кратных
интегралов, когда кубатурные формулы для достижения малой погрешности слишком
громоздки и требуют большого объема в
ы
числений.


П
ри использовании квадратурных или куб
атурных формул, число опер
а
ций быстро
возрастает с ростом раз
мерности интеграла. Например, если для вычисления од
-
номерного интеграла методом трапеций с заданной точн
о
стью необходимо вычислить
сумму порядка
N

слагаемых, то для вычисл
е
ния двойного интеграл
а тем же методом
необходимо сложить порядка
N
2

слагаемых, а для тройно
го интеграла число слагаемых
составляет поря
д
ка
N
3
.

Число испытаний
N
, требующихся для достижения заданной точности
ε

приближе
н
ного значения, в методе

Монте
-
Карло есть величина порядка

и
не
зависит от размерности интеграла
.

Применяется следующий

критерий выбора между кубатурной форму
лой
р
-
го порядка
точности и методом Монте
-
Карло для вычисления с точностью
ε

кратного интеграла
функции
m

пер
е
менных
:

1)

е
сли число

измерений

m

2
р
,
лучше использовать кубатурные или ква
д
ратурные
формулы
;

2)
если

m

� 2
р



ме
тод Монте
-
Карло
.

Например, если
р

 1, тройные интегралы выгоднее вычислять методом Монте
-
Карло,
а одномерные


квад
ратурными формулами.

Если

р

 2, лучше вычис
лять методом Монте
-
Карло пя
тимерные интегр
а
лы, а
одномерные, двойные и трой
ные


квадрату
р
ными или кубатурными формулами.

Рассмотрим конкретные формулы метода Монте
-
Кар
ло для вычисления кратных
интегралов, получающиеся способом, который применялся для в
ы
вода формулы 9.7.

Пусть требуется вычислить двойной инт
е
грал

.


(9.9)

Проведем серию из
N

испытаний случайной точки 
x
i
,
y
i
, где
x
i

равноме
р
но
распределены на отрезке [
a
,
b
],
a

y
i

равномерно ра
с
пред
елены на отрезке [
с
,
d
]. Вычислим
инт
е
грал 9.9 по формуле

=

(9.10)


Для тройного интеграла аналогично пол
у
чим формулу


=

(9.11)


где
x
i

рав
номерно распределены на отрезке [
a
,
b
],
y
i



на отрезке [
с
,
d
],
a

z
i



на отрезке [
р
,
q
];
N



число испы
таний.

Для
m
-
кратного интеграла формула метода Монте
-
Карло имеет вид



.

(9.12)



45. Ре
шение систем линейных уравнений
методом Монте
-
Карло

Рассмотрим применение метода Монте
-
Карло для ре
шения системы л
и
нейных
ура
в
нений


,
i

= 1, ...,
n
.

(9.13)

Пусть система 9.1 приведена к виду

,
i

= 1, ...
,
n
.

(9.14)

где норма матрицы
α

=
||

a
ij

||

удовлетворяет условию ||
α
||  1. Тогда система 9.14 имеет
единственное реш
е
ние.

Приведем без доказательства схематич
е
ское описание метода Монте
-
Карло для
решения си
с
темы линей
ных уравнений 9.14.

Подберем такие множители
v
ij
, чтобы р
е
шения
p
ij

уравнений
α
ij

=
p
ij

v
ij

удовлетворяли
услов
и
ям:

1)

p
ij



0, причем
p
ij

 0, если
α
ij



0;

2)


i

=1, ...,
n
.


Положим

p
i
,

n
+1

=
l



,
i

= 1,

...,


n
;

p

n
+1,
j

= 0,
если

j


n

+
l
;

p
n
+1,
n
+1

=
l
.

Будем трактовать число
p
ij

как вероятность перехода некоторой частицы из состояния
S
i

в
состояние
S
j
. Всего состо
я
ний
п

+

1:

S
1
,
S
2
, ... ,
S
n
+1
,

причем граничное состояние
S
n
+1

=
Г

соответствует ос
тановке частицы, так как
р
п
+
1
,
j

 0, т. е.
частица не мо
жет выйти из с
о
стояния
S
n
+1
.

Матрица с элементами


p
ij


называется
матрицей пере
хода

с
о
стояний {
S
i
}
:

П

. (9.15)

Назовем
траекторией

T
i

совокупность сост
ояний, че
рез которые проходит блуждающая
частица, начиная с состояния
S
i

и заканчивая граничным состоян
и
ем
Г
. Про
межуточные
состо
я
ния частицы обозначим

;

T
i

={
=

Г

}.

(9.16)

Поставим

в соответствие траектории
T
i

случайное чис
ло
X
i
:

X
i


=
β
i

+


+

+




+


, (9.17)

г
де




свободные

члены

системы

(9.14)
с

индексами, совпада
ю
щими с
индексами состояний, обра
зующих траекторию 9.16.

Теорема
9.1
. Математические ожидания случайных величин 9.17 равны корням
системы 9.14:
M
(
X
i
) =
x
i
,
i

= 1, ... ,
п
.


Сформулируем
алг
оритм решения системы линейных уравнений
(9.13)

мет
о
дом
Монте
-
Карло
.

1.

Привести систему 9.1 к виду 9.14:

,
i

= 1, ...,
n
,
||
α
|| 1
.


(9.1
8
)

2.

Подобрать такие числа
v
ij
, чтобы решения
p
ij

урав
нен
ий
α
ij

=
p
ij

v
ij

удовл
е
творяли
услов
и
ям:

1)

p
ij



0, причем
p
ij

 0, если
α
ij



0;

2)


i

=1, ...,
n
.


(9.19)


Вычислить
(
п

+ 1
)
-
й столбец и
(
п

+ 1
)
-
ю строку матрицы перехода:

p
i
,

n
+1

=
l



,
i

= 1, ...,


n
;

p

n
+1,
j

= 0,
если

j


n

+
l
;

p
n
+1,
n
+1

=
l
.


(9.20)



3.

Выполнить для каждого


i

=

1, ... ,
п
:

3.0.


присвоить значение

x
i

= 0;

3.1.


для каждого
k

= 1, ... ,
N

где
N

число испытаний случайной ве
личины
X
i
 выполнить
статистические ис
пытания по алгоритму:

3.1.0.
присвоить

m

= 0;
i
m

=
i
;
v

= 1;
x
i

=
x
i

+
β
i
;

3.1.1.
присвоим переменной очередное значение:
m

=

m

+ 1;

3.1.2.
возьмем случайное число
ξ

из интервала 0; 1, и выберем значение
i
m

по
сл
едующему правилу:

i
m

=
, 1


j



n

3.1.3.

если


j


n

 1, перейти к .1.1;

3.1.4.
для каждого

l


= 0, 1, ... ,
m



2


вычислить

v

=
v

,


x
i

=
x
i

+
v

;

3.2.
вычислить
x
i

=
x
i

/
N
.




Приложенные файлы

  • pdf 17910745
    Размер файла: 2 MB Загрузок: 0

Добавить комментарий