Численные методы решения систем уравнений

1.

Численные методы анализа.
ч.3-4.
«Всё опыт, опыт! Опыт – это вздор.
Значенья духа опыт не покроет.
Всё, что узнали до сих пор,
искать не стоило. И знать не стоит.»
Монолог Бакалавра. «Фауст», Гёте

2.

3. Численные методы решения систем уравнений
3.1. Основные положения
1.
Точные методы – конечные алгоритмы для вычисления корней системы.
2.
Итерационные методы – решение системы путем сходящихся
итерационных процессов.
Источники погрешностей: округления (даже в точных методах) и погрешности
метода.

3.

3.2. Метод Крамера (решение систем линейных уравнений с
помощью обратной матрицы)
a11 x1 a12 x2 ...a1n xn b1
a21 x1 a22 x2 ...a2 n xn b2
...........................................
an1 x1 an 2 x2 ...ann xn bn
Неособенная матрица А:
a11.a12 ...a1n
a21.a22 ...a2 n
A x b
x1
x2
................... ....
an1.an 2 ...ann xn
b1
b2
....
bn
det A , 0
A11. A21... An1
Обратная матрица
А-1:
1 ~ 1 A12 . A22 ... A2 n
A A
...................
1
A1n . A2 n ... Ann
Присоединенная матрица Ã – транспонированная матрица, составленная из
миноров Aij со своими знаками.
x A 1 b
3

4.

3.3. Метод Гаусса (метод гауссовых исключений)
a11 x1 a12 x2 ...a1n xn b1
a21 x1 a22 x2 ...a2 n xn b2
...........................................
an1 x1 an 2 x2 ...ann xn bn
det A 0
a11.a12 ...a1n
x1
b1
a11.a12 ...a1n .b1
a21.a22 ...a2 n
x2
b2
a21.a22 ...a,2 nb2
................... ....
an1.an 2 ...ann xn
....
bn
...................
an1.an 2 ...annbn
x1
x2
....
xn
1
a11.a12 ...a1n .b1
1.. 12 ... 1n .. 1
1..0.....0.. 1
1..0.....0..x1
a21.a22 ...a2 n .b2
0..1..... 2 n .. 2
0..0.....0.. 2
0..0.....0..x2
......................
an1.an 2 ...ann .bn
......................
0..0......1.... n
..................
0..0.....1.. n
..................
0..0.....1..xn
4

5.

3.4. Встроенная функция Lsolve в пакете MathCad
5

6.

3.5. Встроенная функция Find в пакете MathCad
6

7.

3.6. Встроенная функция rref в пакете MathCad
7

8.

3.7. Метод итераций
n
Дана система n линейных уравнений
a x
j 1
ij
j
bi
Предполагается, что диагональные коэффициенты отличны от нуля
aii 0
x1 1 22 x2 ... 1n xn
Система приводится к виду
x2 2 21 x1 23 x... 2 n xn
................................................
xn n n1 x1 n 2 x2 ... n ( n 1) xn 1
или в матричном виде
x β α x
Нулевое приближение
x (0) β
k-ое приближение
x ( k ) β α x ( k 1)
8

9.

Теорема сходимости итерационного ряда
Теорема. Процесс итерации для линейной системы уравнений сходится к
x β α x
единственному ее решению, если какая-нибудь каноническая норма матрицы α
меньше единицы, т.е. достаточным условием сходимости является неравенство
║α ║< 1.
Следствие 1. Процесс итерации сходится, если:
1. Неопределенная норма или m-норма:
2. L1 норма или l-норма:
3. Евклидова норма или k-норма:
l
max
j
m
n
i 1
k
max
ij
i
j 1
ij
1
1
max
i
n
n
i , j 1
2
ij
1
9

10.

Следствие 2. Процесс итерации сходится, если выполнены неравенства :
1.
n
aii ' aij
j 1
2.
n
a jj ' aij
j 1
Условия окончания итерационного процесса:
x ( k 1) x ( k )
x
( k 1)
10

11.

Вычисление норм матриц в пакете MathCad
11

12.

3.8. Метод Зейделя (модификация метода итераций)
При вычислении (k+1) приближения неизвестной xi учитываются уже вычисленные
ранее (k+1) приближения неизвестных x1, x2,… xi-1.
x1 1 22 x2 ... 1n xn
Система приводится к виду
В матричном виде
x2 2 21 x1 23 x... 2 n xn
................................................
xn n n1 x1 n 2 x2 ... n ( n 1) xn 1
x β α x
Нулевое приближение
x
(0)
β
( k 1)
1
x
x2( k 1)
(k+1)-ое приближение
xm( k 1)
xn( k 1)
1 a1 j x
j 2
n
( k 1)
(k )
2 21 x1 a3 j x j
j 3
m 1
n
m amj x (jk 1) amj x (jk )
j 1
j m 1
n 1
n anj x (jk 1)
j 1
n
(k )
j
12

13.

3.9. Метода итераций для систем нелинейных уравнений
Задана система уравнений и начальные приближения
корней x(0)0, x(0)2, …x(0)n:
Система приводится к виду
x1 1 ( x1 , x2 ,...xn )
x2 2 ( x1 , x2 ,...xn )
............................
xn n ( x1 , x2 ,...xn )
x1( k 1)
x1(1) 1 ( x1( 0 ) , x2( 0 ) ,...xn( 0 ) )
(1)
(0)
(0)
(0)
x 2 ( x1 , x2 ,...xn )
1-ое приближение 2
.......................................
(1)
(0)
(0)
(0)
(k )
(k )
(k )
x
(
x
,
x
,...
x
1 ( x1 , x2 ,...xn )
n
n
1
2
n )
2 ( x , x ,...x )
(k+1)-ое приближение x
..........................................
xn( k 1) n ( x1( k ) , x2( k ) ,...xn( k ) )
( k 1)
2
F1 ( x1 , x2 ,...xn ) 0
F2 ( x1 , x2 ,...xn ) 0
............................
Fn ( x1 , x2 ,...xn ) 0
(k )
1
(k )
2
(k )
n
13

14.

4. Интерполирование функций
4.1. Интерполяционная формула Лагранжа
Постановка задачи.
1. На отрезке [a,b] заданы n+1 значения аргумента x : x0, x1,….xn (узлы) и значения
функции yi=f(xi). Требуется построить полином Ln(x) степени не выше n,
имеющий в заданных узлах те же значения, что и функция f(x), т.е. Ln(xi) = yi
при (i=0,1,2,…n).
2. Шаг интерполяции hi+1=(xi+1-xi) может быть произвольным (узлы не являются
равноотстоящими.)
Представим полином в виде
Ln ( x) a0 a1 x ... an x n
В результате получим систему
n+1 уравнений с n+1 неизвестными
a0, a1, a2,…an:
a0 a1 x0 a2 x0 ...an x0 y0
a0 a1 x1 a2 x1 ...an x1 y1
.............................................
a0 a1 xn a2 xn ...an xn yn
14

15.

Неизвестные ai можно найти методом Крамера
Тогда полином примет вид
Qi ( x)
n
Ln ( x) y0Q0 ( x) y1Q1 ( x) ... ynQn ( x)
Функция Qi(x) должна удовлетворять условиям
Её явный вид
an
Qi ( x j ) ij
( x x0 )( x x1 )...( x xi 1 )( x xi 1 )...( x xn )
( xi x0 )( xi x1 )...( xi xi 1 )( xi xi 1 )...( xi xn )
Интерполяционная формула Лагранжа
n
Ln ( x) yn
i 0
( x x0 )( x x1 )...( x xi 1 )( x xi 1 )...( x xn )
( xi x0 )( xi x1 )...( xi xi 1 )( xi xi 1 )...( xi xn )
15

16.

Интерполяционная формула Лагранжа в пакете MathCad
16

17.

4.2. Интерполяционные формулы Ньютона
Рассматриваемые значения аргумента являются равноотстоящими, т.е.
образуют арифметическую прогрессию (шаг интерполяции h = const) .
Определения.
1. Конечные разности первого порядка: Δyi = yi+1 – yi
2. Конечные разности второго порядка: Δ2yi = Δyi+1 – Δyi= yi+2 – 2 yi+1 + yi
3. Конечные разности n-ого порядка: Δnyi = Δn-1yi+1 – Δn-1yi
17

18.

4.2.1. Первая интерполяционная формула Ньютона
Пусть для функции, заданной таблично с постоянным шагом, составлена таблица
конечных разностей. Будем искать интерполяционный полином в виде
Pn ( x) a0 a1 ( x x0 ) a2 ( x x0 )( x x1 )... an ( x x0 )( x x1 )...( x xn 1 )
Коэффициенты a0, a1, ….an найдем из условия совпадения значений функции и
интерполяционного полинома в узлах интерполяции.
Полагая x=x0, находим y0=Pn(x0)= a0.
Далее подставляя x=x1, находим y1=Pn(x1)= a0+a1(x1-x0)= a0+a1h,
x=x2, находим y2=Pn(x2)= a0+a1(x1-x0)+ a2(x2-x0 ) (x2-x1) = a0+2a1h +2a2h2,
18

19.

Найдем коэффициенты a1, ….an :
y0
a1
h
2 y0
a2
2!h 2
m y0
am
m!h 2
y0
2 y0
n y0
Pn ( x) y0
( x x0 )
( x x0 )( x x1 )...
( x x0 )( x x1 )...( x xn 1 )
2
2
h
2!h
n!h
Используя понятие обобщенное степени : x[n]=x(x-h)(x-2h)…[x-(n-1)h] получим
2
y0
n y0
1 y0
2
n
Pn ( x) y0
( x x0 )
(
x
x
)
...
(
x
x
)
0
0
h
2!h 2
n!h 2
Введем переменную q=(x=x0)/h, (q – число шагов).
Тогда первая интерполяционная формула Ньютона примет вид:
Pn ( x) y0 q y0
q (q 1) 2
q (q 1)...(q n 1) n
y0 ...
y0
2!
n!
19

20.

4.2.2. Вторая интерполяционная формула Ньютона
Для интерполирования в конце таблицы применяют вторую интерполяционную
формулу Ньютона :
yn 1
2 yn 2
n y0
Pn ( x) yn
( x xn )
( x xn )( x xn 1 )...
( x xn )...( x x1 )
1!h
2!h 2
n!h n
Введем переменную q=(x-xn)/h, тогда
q (q 1) 2
q (q 1)...(q n 1) n
Pn ( x) yn q yn 1
yn 2 ...
y0
2!
n!
20

21.

4.3. Кубическая сплайн-интерполяция
Кубическим интерполяционным сплайном, соответствующим данной
функции f(x) и данным узлам xi, называется функция S(x),
удовлетворяющая следующим условиям
1. На каждом сегменте [xi-1, xi], i=1,2,…n функция S(x) является
полиномом третьей степени.
2. Функция S(x) , а также ее первая и вторая производные непрерывны на
отрезке [a,b].
3. S(xi)=f(xi), i=0,1,2,…n .
На каждом из отрезков [xi-1, xi], i=1,2,…n будем искать функцию
S(x)= Si(x) в виде полинома третьей степени
21

22.

Si ( x) ai bi ( x xi 1 ) ci ( x xi 1 ) 2 d i ( x xi 1 ) 3 ,
xi 1 x xi
Коэффициенты ai, bi, ci, di подлежат определению (т.е. нахождению) на
всех n элементарных отрезках [xi-1, xi] (i=1,2,…n).
Для того, чтобы система алгебраических уравнений имела решение,
необходимо составить 4n уравнений
Первые 2n уравнений получаются из условия , что график пройдет через
заданные точки
Si ( xi 1 ) ai y i 1
2
3
Si ( xi ) ai bi hi ci hi d i hi
где
hi xi xi 1
i 1,2,...n
22

23.

Следующие (n-1) уравнений вытекают из условия непрерывности первых
производных в узлах интерполяции
d
d
Si 1 ( xi )
Si ( xi )
dx
dx
d
Si ( x) bi 2ci ( x xi ) 3d i ( x xi ) 2
dx
d
Si 1 ( x) bi 1 2ci 1 ( x xi ) 3d i 1 ( x xi ) 2
dx
i 1,2,...n 1
Приравнивая правые части в точка x=xi, получаем
bi 1 bi 2ci hi 3d i hi
2
i 1,2,...n 1
23

24.

Следующие (n-1) уравнений вытекают из условия непрерывности вторых
производных в узлах интерполяции
d2
d2
S ( xi ) 2 Si ( xi )
2 i 1
dx
dx
d
Si ( x) 2ci 6d i ( x xi )
dx
d
Si 1 ( x) 2ci 1 6d i 1 ( x xi )
dx
i 1,2,...n 1
Приравнивая правые части в точка x=xi, получаем
ci 1 ci 3d i hi
i 1,2,...n 1
24

25.

На данном этапе у нас имеется 4n неизвестных и (4n-2) уравнений.
Оставшиеся 2 уравнения можно получить из условия нулевой кривизны
линии в концевых точках (условие свободного закрепления концов). Нулевая
кривизна означает равенство нулю вторых производных в этих точках.
d2
S
(
x
)
0
1
0
dx 2
2
d
S n ( xn ) 0
2
dx
c1 0
2cn 6d n hn 0
Приравнивая правые части в точка x=xi, получаем
ci 1 ci 3d i hi
i 1,2,...n 1
25

26.

В результате получим
ai yi 1
i 1,2,...n
yi yi 1 hi
bi
ci 1 2ci
hi
3
ci 1 ci
di
3hi
i 1,2,...n 1
y yi 1 yi 1 yi 2
hi 1ci 1 2 hi 1 hi ci 3 i
hi 1
hi
c1 0
dn
cn
3d n
i 2,...n
yn yn 1 2hn
bn
cn
hn
3
26

27.

The end
27