Главная страница

Курсовая работа по численным методам для ii курса методические указания СанктПетербург 2009


Скачать 317.74 Kb.
НазваниеКурсовая работа по численным методам для ii курса методические указания СанктПетербург 2009
Дата14.09.2021
Размер317.74 Kb.
Формат файлаpdf
Имя файлаkursovik_2.pdf
ТипКурсовая
#231964

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
Факультет прикладной математики – процессов управления
А. П. ИВАНОВ, И. В. ОЛЕМСКОЙ
КУРСОВАЯ РАБОТА
ПО ЧИСЛЕННЫМ МЕТОДАМ ДЛЯ II КУРСА
Методические указания
Санкт-Петербург
2009

Перейти к оглавлению на странице: 26
КУРСОВАЯ РАБОТА ПО ЧИСЛЕННЫМ МЕТОДАМ
1. Методом Рунге-Кутты IV порядка проинтегрировать диффе- ренциальное уравнение dy dx
+ a · y = ϕ(x),
ϕ(x) = sin kx,
k =

4
,
(1)
на отрезке x ∈ [0; 5] , приняв π = 3.14159265 , с начальным условием y(0) = 0 . Вывести значения интегрируемой функ- ции в точках x i
= 0.5 · i, i = 0, 10.
2. Используя представление искомого решения в виде y(x) =
Z
x
0
e a(t−x)
ϕ(t) dt,
(2)
получить значения y(x) с точностью ε в заданных точках x
i
= 0.5 · i , i = 0, 10 , применив для вычисления интеграла составную квадратурную формулу Симпсона.
3. Найти с заданной точностью ε ближайший к нулю корень уравнения y(x) = 0 .
4. Взяв значения z i
= y(i) + (−1)
i
ε
i
, i = 0, 10 , построить полином наилучшего среднеквадратичного приближения пя- той степени ¯
Q
5
(x) , то есть найти параметры {a i
} полинома
Q
5
(x) = a
0
x
5
+ a
1
x
4
+ a
2
x
3
+ a
3
x
2
+ a
4
x + a
5
из условия
10
X
i=0
( ¯
Q
5
(x i
) − z i
)
2
= min
{a i
}
10
X
i=0
(Q
5
(x i
) − z i
)
2
(3)
Решить эту же задачу, заменив {z i
} точными значениями функции y(x) в точках x i
= i, i = 0, 5 (см. формулу
(1.7)
).
В заданиях 1, 2, 4 построить графики соответствующих функций.
2

Перейти к оглавлению на странице: 26
ГЛАВА 1.
ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗА-
ДАЧИ
КОШИ
ДЛЯ
ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
§1.
Постановка задачи
Пусть в области
D = {a 6 x 6 b, |y i
− y i
0
| 6 b i
} ∈ R
n+1
определена функция f ≡ f (x, y
1
, . . . , y n
),
(x, y) ∈ D.
Необходимо найти решение dy dx
= f (x, y)
(1.1)
удовлетворяющее начальному условию y(x
0
) = y
0
(1.2)
§2.
Одношаговые методы
Одношаговые методы – это методы, которые последователь- но дают приближения y m
к значениям точного решения y(x m
) в каждом узле x m
сетки на основе известного приближения y m−1
к решению в точке x m−1
В общем виде их можно представить :
y m+1
= F (f, x m+1
, x m
, y m+1
, y m
).
(2.1)
Надо отметить, что одношаговые методы вида
(2.1)
в правой ча- сти содержат искомое значение y m+1
. Это позволяет ввести еще один параметр классифицирующий численные методы. Наличие в правой части представления метода в виде
(2.1)
зависимости от искомой функции y m+1
делает его неявным. Отсутствие этой за- висимости – явным, вида:
y m+1
= F (f, x m+1
, x m
, y m
).
(2.2)
3

Перейти к оглавлению на странице: 26
Общая схема явного m – этапного одношагового метода Рунге-
Кутты имеет вид y(x + h) ≈ y(x) +
m
X
i=1
p i
k i
(h),
(2.3)
k i
(h) = hf (x + α
i h, y(x) +
i−1
X
j=1
β
ij k
j
), i = 1, m.
(2.4)
Он идеально приспособлен для практического расчета: во-первых,
не требует вычисления дополнительных начальных значений; во- вторых, позволяет легко менять шаг интегрирования. Функция
Ψ(h) = %
x+h
= y(x + h) − y(x) −
m
X
i=1
p i
k i
(h) ≈ O(h s+1
),
(2.5)
называется методической (локальной) погрешностью одно- шагового метода (методической погрешностью на шаге),
а s – порядком метода.
Выбор постоянных (параметров метода)
p i
, α
i
, β
ij произво- дится так, чтобы разложение методической погрешности
(2.5)
по степеням h начиналось со степени s + 1 при произвольной функ- ции f (x, y).
Расчётная схема четырёхэтапного метода
Рунге-
Кутты четвёртого порядка. (Правило одной шестой).
y(x + h) ≈ y(x) +
1 6
[k
1
(h) + 2k
2
(h) + 2k
3
(h) + k
4
(h)].
(2.6)
где k
1
(h) = hf (x, y(x)),
k
2
(k) = hf (x +
1 2
h, y(x) +
1 2
k
1
(h)),
k
3
(k) = hf (x +
1 2
h, y(x) +
1 2
k
2
(h)),
k
4
(k) = hf (x + h, y(x) + k
3
(h)).
Погрешность расчётной схемы
(2.6)
, как и для всех формул Рунге-
Кутты четвёртого порядка точности, представляется в следующей
4

Перейти к оглавлению на странице: 26
форме:
%
x+h
=
Ψ
(5)
(0)
120
h
5
+ o(h
5
).
(2.7)
Классификация погрешностей.
Требуется найти функцию y(x), которая является решением диф- ференциального уравнения dy dx
= f (x, y(x)),
x ∈ [a, b],
(2.8)
и принимает в точке x
0
= a некоторое определённое значение y(x
0
) = y
0
(2.9)
• В случае, если начальное условие y(x
0
) вычислено (задано) с погрешностью, то вместо задачи
(2.8)
,
(2.9)
, решается задача:
dy
0
dx
= f (x, y
0
(x)),
x ∈ [a, b],
(2.10)
y
0
(x
0
) = ¯
y
0
(2.11)
с изменёнными начальными условиями y
0
− ¯
y
0
= R
0 6= 0.
(2.12)
Решение задачи
(2.10)
,
(2.11)
зависит от ¯
y
0
и не совпадает с решением y(x) исходной задачи
(2.8)
,
(2.9)
Разность
ξ
n
= y(x n
) − y
0
(x n
)
называется неустранимой по- грешностью решения y
0
(x) .

Разность между значением решения y
0
(x n
) задачи
(2.10)
,
(2.11)
и его приближённым значением y
n
,
полученным по формуле y
n
= F (f ; h; x n−1
; y n−1
)
(2.13)
одношагового метода,
ε
n
= y
0
(x n
) − y n
(2.14)
5

Перейти к оглавлению на странице: 26
называется погрешностью метода, или глобальной погрешно- стью метода.
Но вследствие ошибок округления и приближённого вычисле- ния правой части f (x, y) дифференциального уравнения вы- числения значений y n
по формуле
(2.13)
выполняются неточ- но. Фактически найденные значения
ˆ
y n
удовлетворяют не соотношению
(2.13)
, а условию
F (f ; h; x n−1
; ˆ
y n−1
) − ˆ
y n
= δ
n
(2.15)
Невязка δ
n называется погрешностью округления на n -ом шаге.

Разность между точным решением y(x n
)
задачи
(2.8)
,
(2.9)
и приближённым фактически найденным значением ˆ
y n
R
n
= y(x n
) − ˆ
y n
(2.16)
называется полной погрешностью приближённого ре- шения

Величина
η
n
= y n
− ˆ
y n
называется вычислительной по- грешностью.
• Из соотношений
(2)
,
(2.13)
,
(2.16)
,
(2)
следует, что
R
n
= ξ
n
+ ε
n
+ η
n
(2.17)
т.е полная погрешность приближённого решения рав- на сумме неустранимой погрешности, погрешности метода и вычислительной погрешности.
Для полной погрешности справедливо представление
R
n
=
n
X
j=1
(%
j
+ δ
j
) · e
R
xn xj f
y
(ψ,¯
y j
(ψ)) dψ
+ R
0
· e
R
xn x0
f y
(ψ,¯
y
0
(ψ)) dψ
(2.18)
Здесь y j
(x) – интегральная кривая, проходящая через точку
(x j
, ˆ
y j
).
6

Перейти к оглавлению на странице: 26
Из
(2.18)
следует, что полная погрешность приближённого ре- шения задачи Коши
(2.8)
,
(2.9)
в точке x n
равна сумме локальных погрешностей на каждом шаге, взятых с коэффициентами e
R
xn xj f
y
(ψ,¯
y j
(ψ)) dψ
(2.19)
Характер отклонения приближённого решения от точного, т.е. эво- люция полной погрешности, зависит от поведения интегральных кривых уравнения
(2.8)
• Если
∂f
∂y
= 0, т.е. правая часть уравнения не зависит от y ,то полная погрешность равна сумме локальных погрешностей:
R
n
= R
0
+
n
X
j=1
(%
j
+ δ
j
)
(2.20)
• Если
∂f
∂y
> 0, т.е. интегральные кривые расходятся, влияние локальных погрешностей, полученных на предыдущих шагах,
возрастает и полная погрешность больше суммы локальных погрешностей.
• Если
∂f
∂y
< 0, т.е. интегральные кривые сближаются, влияние локальных погрешностей, полученных на предыдущих шагах,
ослабевает и полная погрешность меньше суммы локальных погрешностей.
Такое поведение характерно как для полной погрешности, так и для отдельных её частей.
7

Перейти к оглавлению на странице: 26
Практическая реализация явных одношаговых мето- дов типа Рунге-Кутты решения задачи Коши.
Полагаем, что в нашем вычислительном процессе погрешно- стью округления и неустранимой погрешностью можно пренебречь.
Оценки погрешностей необходимы, с одной стороны, чтобы обеспечить длину шага h, достаточно малую для достижения тре- буемой точности вычисляемых результатов, а с другой стороны –
чтобы гарантировать достаточно большую длину шага во избежа- ние бесполезной вычислительной работы.
Метод Рунге оценки полной погрешности.
Полагаем, что в точке x n
по m -этапному методу (
2.3
) s -ого порядка точности с постоянным шагом h вычислено приближён- ное решение ¯
y n
исходной задачи Коши. Считаем, что для полной погрешности метода Рунге-Кутты s− о порядка справедливо пред- ставление:
y(x n
) − ¯
y n
≈ z(x n
)h s
,
(2.21)
где z(x n
)h s
− главный член аcимптотического разложения полной погрешности. Используя ту же расчётную формулу с шагом h
2
,
вычислим в той же точке x n
другое значение решения ˜
y n
. Для этого потребуется в два раза больше шагов и погрешность прибли- жённого решения в этом случае может быть представлена y(x n
) − ˜
y n
≈ z(x n
)
 h
2

s
(2.22)
Исключим из
(2.21)
,
(2.22)
точное значение решения y(x n
) И пра- вило Рунге оценки глобальной погрешности может быть записано следующим образом:
¯
R
n
= y(x n
) − ¯
y n


y n
− ¯
y n
)
(1 − 2
−s
)
(2.23)
˜
R
n
= y(x n
) − ˜
y n


y n
− ¯
y n
)
(2
s
− 1)
(2.24)
В качестве решения в точке x n
примем значение
˜
y n
как более точное по сравнению с
¯
y n
. Для него имеем оценку погрешности
8

Перейти к оглавлению на странице: 26
(2.24)
. Эта величина может быть как больше, так и меньше некото- рого значения ε, являющегося наперёд заданной допустимой погрешностью.
Если | ˜
R
n
| 6 ε, то заданная точность приближённого решения достигается. Полученное приближённое решение можно уточнить,
прибавив к нему величину главного члена погрешности, т.е. поло- жив y(x n
) ≈ y n
= ˜
y n
+ ˜
R
n
При этом y(x n
) − y n
= O(h s+1
).
9

Перейти к оглавлению на странице: 26
ГЛАВА 2.
ВЫЧИСЛЕНИЕ ОПРЕДЕЛЕННОГО ИНТЕ-
ГРАЛА
§1.
Квадратурные формулы
Будем рассматривать задачу о вычислении однократного ин- теграла
J (F ) =
Z
b a
F (x) dx
(1.1)
с помощью конечного числа значений интегрируемой функции.
Интервал интегрирования
[a, b]
есть любой конечный (или бесконечный) отрезок числовой оси. Подынтегральная функция
F (x) – любая интегрируемая в смысле Римана функция.
Каждое правило основано на замене интегрируемой функции на какую-либо элементарную функцию – алгебраический много- член, рациональную функцию, тригонометрический многочлен и др. Чтобы такая замена давала хорошую точность, требуется, что- бы заменяемая функция F (x) обладала высоким порядком глад- кости.
Если F (x) имеет какие-нибудь особенности, мы заинтересо- ваны в выделении их. Выделение делается обычно при помощи разложения F (x) на два сомножителя F (x) = p(x) · f (x), где
• p(x) имеет особенности того же типа, что и F (x), и назы- вается весовой функцией
;
• f (x) – гладкая функция.
Такое разложение приводит
(1.1)
к виду:
J (F ) =
Z
b a
p(x)f (x) dx.
(1.2)
Считая вес p(x) фиксированным, а f (x) любой гладкой функцией,
будем строить правило интегрирования, рассчитанное на функции,
имеющие одинаковые, заранее известные особенности. С учётом
10

Перейти к оглавлению на странице: 26
всего вышесказанного будем строить правила вычисления опреде- лённого интеграла следующего вида:
J (F ) =
Z
b a
p(x)f (x) dx ≈
n
X
k=1
A
k f (x k
) = S
n
,
x k
∈ [a, b],
(1.3)
где
• f (x) ∈ Φ(a, b) – некоторый класс функций f (x) , определён- ных на [a, b] ;
• p(x) –
весовая функция
– некоторая фиксированная неот- рицательная на [a, b] функция, для которой
Z
b a
p(x) dx > 0,
(1.4)
и для ∀f (x) ∈ R существует
Z
b a
p(x)|f (x)| dx.
(1.5)
Определение 1.1. Равенство
(1.3)
называют формулой механи- ческих квадратур или просто квадратурной формулой, а S
n
=
P
n k=1
A
k f (x k
) – квадратурной суммой; A
k
– квадратурными коэф- фициентами формулы
(1.3)
; x k
– узлами квадратурной формулы
(1.3)
Будем предполагать, что узлы квадратурной формулы
(1.3)
упорядочены по возрастанию a = x
0
< x
1
< . . . < x n
=
b , n, x k
, A
k
,
k = 1, . . . , n; – являются параметрами квадратур- ной формулы
(1.3)
и их следует выбирать так, чтобы достигнуть
“возможно лучшего” результата интегрирования для всех функций избранного класса.
Определение 1.2. Величина
R
n
(f ) =
Z
b a
p(x)f (x) dx −
n
X
k=1
A
k f (x k
) = J (f ) − S
n
(f ).
(1.6)
называется методической погрешностью (остаточным членом)
квадратурной формулы
(1.3)
11

Перейти к оглавлению на странице: 26
При зтом R
n
(f ) зависит от свойств функции f и от выбора квадратурной формулы, т. е. от узлов и коэффициентов.
Определение 1.3. Квадратурная формула
(1.3)
имеет алгебраиче- скую степень точности m, если она верна для любых многочленов степени m и не верна для многочленов степени m + 1 .
Это равносильно тому, что равенство
Z
b a
p(x)x i
dx =
n
X
k=1
A
k x
i k
,
(1.7)
выполняется для i = 0, 1, . . . , m и не выполняетсядля i = m + 1.
Или, что то же самое,
R
n
(x i
) = 0,
i = 0, 1, . . . , m;
R
n
(x m+1
) 6= 0.
(1.8)
§2.
Квадратурная формула Симпсона
Трёхточечная квадратурная формула – формула Симпсона –
имеет вид :
Z
b a
f (x) dx ≈
(b − a)
6

f (a) + 4 · f (
a + b
2
) + f (b)

(2.1)
Она является точной для всех многочленов третьей степени, т.е. ее алгебраическая степень точности равна трем.
В классе C
4
(a,b)
остаточный член квадратурной формулы
(2.1)
R
3
(f ) =
Z
b a
f
(4)
(ξ(x))
4!
(x − a)(x −
a + b
2
)
2
(x − b) dx =
= −
1 90
(
b − a
2
)
5
f
(4)
(η),
η ∈ [a, b].
(2.2)
Составные квадратурные формулы.
Основная идея мето- да заключается в том, что для повышения точности интегрирова- ния отрезок [a, b] делят на несколько частей. Применяют избран- ную квадратурную формулу к каждой отдельной части и резуль- таты складывают.
12

Перейти к оглавлению на странице: 26
Для большинства квадратурных формул методическая по- грешность
R
n
(f )
зависит от величины отрезка интегрирования и может быть представлена в виде:
R
n
(f ) = (b − a)
k
C(a, b),
(2.3)
где C(a, b) есть медленно меняющаяся функция от a,
b и k.
Эта зависимость показывает, что если уменьшить отрезок интегри- рования в p раз, то R
n
(f ) уменьшится в p k
раз. Чем больше k, тем значительнее будет уменьшение.
Разобъём исходный интервал интегрирования на частичные отрезки:
[a, b] : a = z
0
< z
1
< z
2
< . . . < z k
= b.
В самом общем случае z i+1
− z i
= h i
> 0 и h i
предполагаются различными.
Теперь применим на каждом из частичных отрезков [z i
, z i+1
]
квадратурное правило (в общем случае своё) для вычисления ин- теграла:
J
(i)
(f ) =
Z
z i+1
z i
p(x)f (x) dx =
n i
X
j=1
A
ij f (x ij
) + R
n i
(f ) =
= S
i n
i
(f ) + R
n i
(f ),
i = 0, k − 1.
(2.4)
Просуммируем по i левые и правые части
(2.4)
:
J (f ) =
Z
b a
p(x)f (x) dx =
k−1
X
i=0
J
(i)
=
k−1
X
i=0
Z
z i+1
z i
p(x)f (x) dx ≈

k−1
X
i=0
n i
X
j=1
A
ij f (x ij
) =
k−1
X
i=0
S
i n
i
(f ) = S
N
(f ),
(2.5)
N =
k−1
X
i=0
n i
;
R
N
= J (f ) − S
N
(f ) =
k−1
X
i=0
R
n i
(f ),
R
n i
(f ) = J
i
(f ) − S
i n
i
(f ).
В качестве квадратурных формул
(2.4)
могут быть использованы квадратурные формулы всех существующих типов. Причём как в
13

Перейти к оглавлению на странице: 26
смешанном составе – в составной квадратурной формуле использу- ется несколько типов, так и в однородном – в составной квадратур- ной формуле используются только квадратурные формулы одного типа. Тоже самое можно сказать и о числе узлов квадратурной фор- мулы используемом на каждом из частичных отрезков - их число может зависеть от номера отрезка. Такой подход увеличения точ- ности применим к простейшим формулам Ньютона-Котеса (здесь
– Симпсона).
Составная квадратурная формула Симпсона.
Разобъём интервал интегрирования [a, b] на k = 2m частичных отрезков,
z i
= a + iH, i = 0, k
H =
b−a k
и будем применять квадратурное правило на сдвоенном частичном отрезке [a + (2j − 2)H, a + (2j)H]
j = 1, m.
В качестве квадратурного правила используемого на каждом из сдвоенных частичных отрезков будет использоваться квадратурное правило Симпсона
(2.1)
:
J
(j)
(f ) =
Z
z
2j z
2j−2
f (x) dx ≈
z
2j
− z
2j−2 6
[f (z
2j−2
) + 4f (z
2j−1
)+
+f (z
2j
)] =
H
3
[f (z
2j−2
+ 4f (z
2j−1
) + f (z
2j
)] ,
j = 1, m. (2.6)
Просуммировав по всем сдвоенным частичным отрезкам
[a, a +
2H],
[a + 2H, a + 4H], . . . выпишем составную квадратурную фор- мулу парабол
Z
b a
f (x) dx =
b − a
3k
[f (z
0
) + f (z k
) + 2(f (z
2
) + f (z
4
) + . . .
+f (z k−2
)) + 4(f (z
1
) + f (z
3
+ . . . + f (z k−1
))] + R
k
(f ).(2.7)
Воспользуемся результатами оценки методической погрешности для квадратурной формулы парабол
(2.2)
:
R
k+1
(f ) = −
H
4 180
Z
b a
f
(4)
(x) dx + o(H
4
) =
(2.8)
= C
4
H
4
+ o(H
4
).;
o(H
4
)
H
4
→ 0, H → 0,
где C
4
≡ C
4
(f ) ≡ −
1 180
Z
b a
f
(4)
(x) dx.
14

Перейти к оглавлению на странице: 26
Практические способы оценки погрешности составных квадратурных формул.
Оказывается, что подобно представле- нию
(2.3)
для погрешности составной квадратурной формулы
(2.5)
для широкого класса функций f (·) можно написать:
R
k
(f ) = C
m
H
m
+ o(H
m
).
(2.9)
где m – натуральное число, C
m
= C
m
(f ) – некоторая констан- та, зависящая лишь от f (·)
и типа формулы, но не зависящая от H,
o(H
4
)
H
4
→ 0, H → 0. Формула
(2.9)
дает ассимптотиче- ское представление (разложение) погрешности
(2.5)
по параметру
H - шагу равномерного разбиения интервала интегрирования. Она справедлива, в частности, когда малые квадратурные формулы
(2.1)
(Симпсона, например,) имеют алгебраическую степень точно- сти m − 1, а функция f (·) – непрерывную (интегрируемую) на
[a, b] производную f m
(·).
Практическая оценка константы C
m проводится следующим образом. Фиксируем два значения шага разбиения –
H
1
= H
и
H
2
= H/L, L > 1,
и проводим расчёты на двух равномерных сет- ках с шагом H
1
и H
2
соответственно по исследуемой квадратур- ной формуле. Предполагаем, что для её методической погрешности верно аcимптотическое разложение
(2.9)
, т.е.
R
(H
1
)
k
(f ) = J (f ) − S
(H
1
)
k
= C
m
H
m
+ o(H
m
)
(2.10)
и
R
(H
2
)
k
(f ) = J (f ) − S
(H
2
)
k
= C
m
 H
L

m
+ o
 H
L

m
(2.11)
Исключая из
(2.10)
,
(2.11)
неизвестное значение интеграла J, най- дём
S
(H
2
)
k
− S
(H
1
)
k
= C
m

H
m

 H
L

m

+ o(H
m
)
(2.12)
а отсюда
C
m
=
S
(H
2
)
k
− S
(H
1
)
k
H
m
(1 − L
−m
)
+
o(H
m
)
H
m
(2.13)
15

Перейти к оглавлению на странице: 26
Подставляя это выражение для C
m в равенства
(2.10)
,
(2.11)
и от- брасывая бесконечно малые o(H
m
), найдём приближённые пред- ставления погрешностей
R
(H
1
)
k
, R
(H
2
)
k составной квадратурной формулы на двух равномерных сетках :
R
(H
1
)
k
(f ) = J (f ) − S
(H
1
)
k

S
(H
2
)
k
− S
(H
1
)
k
1 − L
−m
(2.14)
R
(H
2
)
k
(f ) = J (f ) − S
(H
2
)
k

S
(H
2
)
k
− S
(H
1
)
k
L
m
− 1
(2.15)
Этот способ оценивания методической погрешности называется правилом Рунге.
Так при
L = 2,
т.е. шаг дробления умень- шается в два раза, погрешность составной квадратурной формулы
R
(H
2
)
k
(f ) с шагом H
2
меньше погрешности R
(H
1
)
k
(f ) составной квадратурной формулы с шагом H
1
почти в 2
m раз.
Используя правило Рунге можно решить вопрос о нахожде- нии оптимального шага разбиения H
o интервала интегрирования,
обеспечивающего (в первом приближении) вычисление интеграла с заданной точностью ε.
ε = |R
(H
o
)
k
(f )| ≈ |C
m
H
m o
| =
S
(H
2
)
k
− S
(H
1
)
k
H
m
(1 − L
−m
)
H
m o
,
(2.16)
а отсюда и получаем порядок оптимального шага разбиения
H
o
≈ H


ε(1 − L
−m
)
S
(H
2
)
k
− S
(H
1
)
k


1
m
(2.17)
16

Перейти к оглавлению на странице: 26
ГЛАВА 3.
МЕТОД НЬЮТОНА РЕШЕНИЯ УРАВНЕ-
НИЙ
§1.
Основы метода
Рассмотрим вопрос о решении скалярного уравнения f (x) = 0,
x ∈ [α, β]
(1.1)
с использованием метода Ньютона в предположении, что решение
¯
x уравнения существует и ¯
x ∈ [α, β] . При условии, что функция f (·) дифференцируема на рассматриваемом интервале и при на- личии хорошего приближения x k
к корню ¯
x функции можно ис- пользовать метод Ньютона, называемый также методом линеари- зации или методом касательных. Его расчётные формулы могут быть получены заменой исходного уравнения
(1.1)
в окрестности корня линейным уравнением (здесь x k
– приближение корня ¯
x )
f (x k
) + f
0
(x k
)(x − x k
) = 0.
(1.2)
Решение этого уравнения принимается за очередное приближение x
k+1
к искомому корню ¯
x исходного уравнения
(1.1)
:
x k+1
= x k

f (x k
)
f
0
(x k
)
(1.3)
Метод Ньютона имеет простую геометрическую интерпретацию:
график функции заменяется касательной к нему в точке (x k
, f (x k
))
и за очередное приближение x k+1
принимается абсцисса точки пе- ресечения этой касательной с осью Ox . Используя эту интерпре- тацию, легко получить расчётные формулы
(1.3)
метода Ньютона и вследствие этой интерпретации он именуется, как указано выше,
методом касательных.
Ясно, что сходимость последовательности {x k
} к корню зави- сит от свойств функции f (·) и не всегда имеет место. Пусть выбра- но начальное приближение x
0
. Легко представить, что уже прибли- жение x
1
не попадает на исходный интервал определения функции и процесс итераций останавливается.
Приведём полезную теорему, гарантирующую сходимость ме- тода Ньютона.
17

Перейти к оглавлению на странице: 26
Теорема 1.1. Если для функции f (·) ∈ C
2
[α, β] выполнены усло- вия
• f (α) · f (β) < 0 ,
• f
0
(x) 6= 0 и f
00
(x) 6= 0 на [α, β] (и, следовательно, сохраняют определённые знаки при x ∈ [α, β]) ,
то исходя из начального приближения x
0
∈ [α, β] , для которого f (x
0
) · f
00
(x
0
) > 0 , по формуле
(1.3)
можно вычислить единствен- ный корень ¯
x уравнения
(1.1)
с любой степенью точности. ♣
Замечание.
Практическим критерием окончания вычислений является выполнение условия |x n+1
− x n
| < ε , где ε – требуемая точность вычисления корня.

Если же условия теоремы
1.1
не выполняются или проверка их затруднительна, то очередное “приближение” x k+1
может оказать- ся вне интервала, на котором расположен корень ¯
x . В этом случае x
k+1
строится либо методом половинного деления либо методом хорд. В первом случае полагают x
k+1
=
a k
+ b k
2
,
(1.4)
во втором –
x k+1
= a k

b k
− a k
f (b k
) − f (a k
)
· f (a k
).
(1.5)
Здесь a k
, b k
– левый и правый конец интервала, которому принад- лежит корень ¯
x на предыдущем шаге.
На начальном этапе полагаем a
0
= α, b
0
= β . Пусть для определённости f (a) < 0, f (b) > 0. Если x
1
∈ [α, β] , то вычислив c = f (x
1
) , полагаем a
1
= c , b
1
= b
0
при c < 0 , и a
1
= a
0
, b
1
= c при c > 0 и повторяем вычисления.
Если же приближение x
1 6∈ [α, β] , то применяем формулы
(1.4)
либо
(1.5)
и поступаем как и выше: вычисляя c = f (x
1
) , полагаем a
1
= c , b
1
= b
0
при c < 0 , и a
1
= a
0
, b
1
= c при c > 0 и применяем метод Ньютона.
18

Перейти к оглавлению на странице: 26
Комментарии к выполнению задания
Формула
(1.3)
для функции, заданной уравнением
(1)
, запи- шется в виде x
k+1
= x k

y(x k
)
y
0
(x k
)
= x k

y(x k
)
ϕ(x k
) − ay(x k
)
(1.6)
Следовательно, для её использования достаточно проинтегриро- вать уравнение
(1)
на отрезке [0, x k
] , либо вычислить интеграл
(2)
в пределах от 0 до x k
. Можно также использовать аналитическое представление для решения уравнения
(1)
, которое имеет вид y(x) =
a · sin kx − k · cos kx + ke
−ax a
2
+ k
2
(1.7)
Начальное приближение для использования формулы
(1.6)
может быть найдено, например, графическим методом, т.к. к моменту ре- шения этой задачи уравнение
(1)
уже проинтегрировано.

19

Перейти к оглавлению на странице: 26
ГЛАВА 4.
ЗАДАЧА АППРОКСИМАЦИИ
Для решения поставленной в пункте 4 задачи использу- ется метод наименьших квадратов. он состоит в следующем.
§1.
Линейная задача метода наименьших квадратов
Пусть на отрезке [a, b] задана некоторая функция f (·) , причём её значения {y i
= f (x i
)} в узлах сетки {x
0
, x
1
, . . . , x n
} известны с погрешностями {ε
i
} , то есть вместо набора значений {y i
} имеем набор {˜
y i
= y i

i
} . (Далее под {y i
} будем понимать заданные зна- чения функции, т.е. с погрешностями, вводя обозначение ˜
y i
лишь в случае необходимости). Пусть, кроме того, на [a, b] определены функции ϕ
j
(·) ∈ Φ , j = 0, m .
Введём в рассмотрение обобщённый полином
P
m
(x) =
m
X
j=0
a j
ϕ
j
(x).
Пусть a – вектор коэффициентов полинома P
m
(·) , y – вектор зна- чений функции f (·) и, наконец, ϕ(·) – вектор-функция, составлен- ная из значений {ϕ
i
(x)} :
a = (a
0
, a
1
, . . . , a m
)
T
,
y = (y
0
, y
1
, . . . , y n
)
T
,
ϕ(x) = (ϕ
0
(x), ϕ
1
(x), . . . , ϕ
m
(x))
T
,
Введём также функции
σ(a, y) =
n
X
i=0
(P
m
(x i
) − y i
)
2
,
δ(a, y) =
r
σ(a, y)
n + 1
Функцию δ(a, y) назовём среднеквадратичным уклонением обоб- щённого полинома P
m
(·) от функции f (·) на системе узлов {x i
} .
Отказываясь теперь от условия P
m
(x i
) = y i
, i = 0, n , поставим задачу так: найти обобщённый полином ¯
P
m
(·) = ¯
a
T
ϕ(·) , для кото- рого среднеквадратичное уклонение минимально :
δ(¯
a, y) = min a
δ(a, y).
20

Перейти к оглавлению на странице: 26
Поставленную здесь задачу называют линейной задачей метода наименьших квадратов или просто методом наименьших квадра- тов (МНК). Если искомый полином существует, то будем называть его многочленом наилучшего среднеквадратичного приближения.
Отметим, что минимум функций σ(·) и δ(·) достигается на одном и том же векторе ¯
a , поэтому фактически будем вести ми- нимизацию функции σ(a, y) , а не δ(·) .
Простейший подход к решению задачи – использование необ- ходимых условий в задаче поиска экстремума для дифференциру- емой функции:
∂σ(a, y)
∂a k
a=¯
a
= 0,
k = 0, m.
(1.1)
В дополнение к уже введённым ранее обозначениям примем также следующее:
Q =




ϕ
0
(x
0
),
ϕ
1
(x
0
),
ϕ
m
(x
0
),
ϕ
0
(x
1
),
ϕ
1
(x
1
),
ϕ
m
(x
1
),
ϕ
0
(x n
),
ϕ
1
(x n
),
ϕ
m
(x n
)




Кроме того будем считать, что под скалярным произведением двух векторов u, v ∈ R
k понимается число
(u, v) = u
T
v = u · v =
k
X
i=1
u i
v i
,
а под нормой векторов – евклидова норма kyk
2
= (y, y) . Тогда в новых обозначениях
σ(a, y) = kQa − yk
2
= (Qa − y, Qa − y) =
= (Qa, Qa) − 2(Qa, y) + (y, y) =
= (a, Q
T
Qa) − 2(a, Q
T
y) + kyk
2
Для определения параметров искомого полинома в соответствии с формулой
(1.1)
имеем СЛАУ:
Ha = b,
где H = Q
T
Q, b = Q
T
y.
(1.2)
21

Перейти к оглавлению на странице: 26
В предложенной в пункте 4 задаче её параметры соотносятся с ис- пользованными выше параметрами следующим образом:
a = (a
0
, a
1
, . . . , a
5
)
T
, y = z = (z
0
, z
1
, . . . , z
10
)
T
, ϕ(x) = (1, x, . . . , x
5
)
T
,
σ(a, y) = σ(a, z) =
10
X
i=0
(Q
5
(x i
) − z i
)
2
, δ(a, z) =
r
σ(a, z)
5
,
Q =




1,
x
0
,
x
5 0
,
1,
x
1
,
x
5 1
,
1,
x
5
,
x
5 5




Остаётся для определения параметров полинома ¯
Q
5
(x) решить си- стему линейных алгебраических уравнений
(1.2)
. В предложенном варианте матрица H системы
(1.2)
неособая и задача имеет един- ственное решение.
Выполняя вторую часть задания из п.4, Вы должны получить так называемый интерполяционный полином ˜
Q
5
(x) , обладающим свойством: ˜
Q
5
(i) = y(i) , i = 0, 5 .
(Построив полином – убедиться!)
§2.
Методы Гаусса
Метод Гаусса относится к точным методам решения систем линейных алгебраических уравнений (СЛАУ). Под точными мето- дами решения СЛАУ понимают методы, которые позволяют за ко- нечное число арифметических операций при отсутствии округле- ний получить точное решение СЛАУ при точно заданной правой части СЛАУ и и точно заданных элементах её матрицы. Перечис- лим здесь некоторые из них с указанием алгоритмов их использо- вания.
Методы Гаусса (или методы исключения неизвестных) опти- мальны для решения СЛАУ общего вида по количеству арифмети- ческих операций, необходимых для нахождения решения этой си- стемы. Пусть СЛАУ записана в виде
Ax = b.
(2.1)
Простой метод Гаусса состоит в следующем: в СЛАУ
(2.1)
, запи- санной в координатной форме, производится деление первого урав- нения системы на коэффициент при первом неизвестном, после че-
22

Перейти к оглавлению на странице: 26
го это преобразованное уравнение, будучи умноженным последо- вательно на все коэффициенты при первом неизвестном из всех нижележащих уравнениях системы, вычитается из соответствую- щих уравнений. Тем самым первое неизвестное оказывается “ис- ключенным” из всех уравнений системы, кроме первого. Оставив на время в покое первое уравнение, таким же образом исключа- ем из системы уравнений с номерами 2–n второе неизвестное и так далее. Результатом проделанной работы явится уравнение, содер- жащее лишь последнее неизвестное системы
(2.1)
, из которого оно и находится. Этот этап преобразований СЛАУ носит название прямо- го хода метода Гаусса. Вслед за найденным неизвестным находятся поочередно и остальные. Этот этап называется обратным ходом ме- тода Гаусса.
Отметим очевидный недостаток простого метода Гаусса – воз- можное обращение в нуль ведущих элементов, т.е. тех элементов,
на которые приходится делить в прямом ходе метода Гаусса.
Легко указать и способ избавления от этого недостатка простого метода
Гаусса – перестановка уравнений системы, которая, очевидно, не меняет решения СЛАУ. Ввиду того, что решается СЛАУ с неособой матрицей, при исключении первого неизвестного в первом столбце матрицы A найдётся ненулевой элемент и строка, содержащая этот элемент, переставляется на место первого уравнения. С целью ми- нимизации погрешности округления обычно выбирают не просто ненулевой элемент в столбце, а элемент, максимальный по модулю и этот процесс повторяется на каждом этапе прямого хода метода
Гаусса. Данная модификация называется методом Гаусса с выбором максимального элемента по столбцу. Можно обобщить этот подход и выбирать максимальный элемент по всей матрице, что приводит к необходимости переобозначения неизвестных и потому использу- ется значительно реже.
Отметим здесь, что если матрица A системы
(2.1)
обладает свойством диагонального преобладания (т.е. для всех строк мат- рицы A = {a ij
} выполняются неравенства |a ii
| >
P
j6=i
|a ii
| ), то обращение в нуль ведущих элементов в методе Гаусса не происхо- дит.
В подробной записи метод Гаусса выглядит так: имеется
23

Перейти к оглавлению на странице: 26
СЛАУ, записанная в виде











a
11
x
1
+ a
12
x
2
+ a
13
x
3
+ . . . + a
1n x
n
= b
1
a
21
x
1
+ a
22
x
2
+ a
23
x
3
+ . . . + a
2n x
n
= b
2
a
31
x
1
+ a
32
x
2
+ a
33
x
3
+ . . . + a
3n x
n
= b
3
a n1
x
1
+ a n2
x
2
+ a n3
x
3
+ . . . + a nn x
n
= b n
(2.2)
Пусть a
11 6= 0, иначе меняем порядок уравнений с тем, что- бы было выполнено данное условие. Разделим первое уравнение на a
11
, после чего оно примет вид:
x
1
+ a
(1)
12
x
2
+ a
(1)
13
x
3
+ . . . + a
(1)
1n x
n
= b
(1)
1
(2.3)
где a
(1)
1j
= a
1j
/a
11
,
b
(1)
1
= b
1
/a
11
Умножим уравнение
(2.3)
на a
21
и вычтем полученное уравне- ние из второго уравнения системы
(2.2)
. Аналогично преобразуем остальные уравнения. В итоге получим СЛАУ









x
1
+ a
(1)
12
x
2
+ a
(1)
13
x
3
+ . . . + a
(1)
1n x
n
= b
(1)
1
a
(1)
22
x
2
+ a
(1)
23
x
3
+ . . . + a
(1)
2n x
n
= b
(1)
2
a
(1)
n2
x
2
+ a
(1)
n3
x
3
+ . . . + a
(1)
nn x
n
= b
(1)
n
,
(2.4)
где введены обозначения a
(1)
2j
=
a
2j
− a
(1)
1j a
21
,
a
(1)
3j
= a
3j
− a
(1)
1j a
31
, . . . ,
a
(1)
nj
=
a nj
− a
(1)
1j a
n1
,
b
(1)
j
= b j
− b
(1)
j a
j1
, . . . ,
j = 2, n.
Теперь, оставив без изменений первое уравнение системы
(2.4)
,
можно применить описанную процедуру к системе из (n − 1) -го уравнений, исключив неизвестное x
2
из третьего и последующих уравнений. Получим систему вида









x
1
+ a
(1)
12
x
2
+ a
(1)
13
x
3
+ . . . + a
(1)
1n x
n
= b
(1)
1
x
2
+ a
(2)
23
x
3
+ . . . + a
(2)
2n x
n
= b
(2)
2
a
(2)
n3
x
3
+ . . . + a
(2)
nn x
n
= b
(2)
n
,
24

Перейти к оглавлению на странице: 26
где a
(2)
2j
= a
(1)
2j
/a
(1)
22
,
b
(2)
2
= b
(1)
2
/a
(1)
22
,
a
(2)
3j
= a
(1)
3j
− a
(2)
2j a
(1)
32
,
a
(2)
nj
= a nj
(1) − a
(2)
2j a
(1)
n2
,
b
(2)
j
= b j
(1) − b
(2)
j a
(1)
j
, . . . ,
j = 3, n.
Продолжая аналогичные вычисления, в итоге получим эквивалент- ную исходной треугольную систему













x
1
+ a
(1)
12
x
2
+ a
(1)
13
x
3
+ . . . + a
(1)
1n x
n
= b
(1)
1
x
2
+ a
(2)
23
x
3
+ . . . + a
(2)
2n x
n
= b
(2)
2
x
3
+ . . . + a
(3)
3n x
n
= b
(3)
3
x n
= b
(n)
n
,
откуда и находятся последовательно все компоненты решения.

25

ОГЛАВЛЕНИЕ
Курсовая работа по численным методам . . . . . . . . . . . . . . .
2
Глава 1. Численные методы решения задачи Коши для обыкновенных дифференциальных урав- нений . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
§ 1. Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
§ 2. Одношаговые методы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Глава 2. Вычисление определенного интеграла . . . . . . .
10
§ 1. Квадратурные формулы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
§ 2. Квадратурная формула Симпсона . . . . . . . . . . . . . . . . . . . . . .
12
Глава 3. Метод Ньютона решения уравнений . . . . . . . . .
17
§ 1. Основы метода . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
Глава 4. Задача аппроксимации . . . . . . . . . . . . . . . . . . . . . . . .
20
§ 1. Линейная задача метода наименьших квадратов . . . . . . . .
20
§ 2. Методы Гаусса . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22 26


написать администратору сайта