Mathcad - ЛР2
.pdfТеперь найдём решение СЛУ методом Якоби:
æ31 |
28 |
19 ö |
æ-12 ö |
|
A = ç19 |
31 |
21 ÷ |
B = ç 21 |
÷ |
ç |
|
÷ |
ç |
÷ |
è75 |
38 |
48 ø |
è 56 |
ø |
check_diag_dom(A) = "Нет диагонального преобладания..."
В данном случае диагональное преобладание отсутствует, следовательно, не выполнено достаточное условие сходимости. Добъемся диагонального преобладания путем эквивалентного преобразования исходной матрицы:
|
æ75 |
38 |
48 ö |
|
æ 56 ö |
|
|
æ75 |
38 |
48 ö |
|
æ 56 ö |
|||
A := |
ç |
19 |
31 |
21 ÷ |
B := |
ç |
21 ÷ |
|
A := |
ç19 |
31 |
21 ÷ |
B := |
ç 21 |
÷ |
|
ç |
|
|
÷ |
|
ç |
÷ |
|
|
ç |
|
÷ |
|
ç |
÷ |
|
è31 28 19 ø |
|
è |
-12 ø |
|
|
è12 |
-3 -2 ø |
|
è-33 |
ø |
||||
|
æ12 -3 -2 ö |
|
|
æ-33 ö |
|
æ12 |
-3 -2 ö |
|
æ-33 ö |
||||||
A := |
ç |
19 |
31 |
21 ÷ |
B := |
ç 21 |
÷ |
A := |
ç 7 |
34 |
23 ÷ |
B := |
ç 54 |
÷ |
|
|
ç |
|
|
÷ |
|
|
ç |
÷ |
|
ç |
|
÷ |
|
ç |
÷ |
|
è75 |
38 |
48 ø |
|
|
è 56 |
ø |
|
è75 |
38 |
48 ø |
|
è 56 |
ø |
|
|
æ12 |
-3 -2 ö |
|
|
æ-33 ö |
|
|
|
|
|
|
|
|||
A := |
ç |
7 |
34 |
23 ÷ |
B := |
ç 54 |
÷ |
|
|
|
|
|
|
|
|
|
ç |
|
|
÷ |
|
|
ç |
÷ |
|
|
|
|
|
|
|
|
è 3 |
56 |
60 ø |
|
|
è254 |
ø |
|
|
|
|
|
|
|
|
check_diag_dom(A) = "Есть диагональное преобладание!"
Теперь, когда диагональное преобладание достигнуто, можно применять итерационные методы Якоби и Зейделя:
e := 10− 6
X_Jacobi := JacobiPlus(A , B , 0 , e)
X - двумерный массив, где каждый из его элементов Xk- вектор-столбец, содержащий значения переменных системы, вычисленных на k-ом шаге.
k := JacobiPlus(A , B , 1 , e) |
|
Количество итераций k = 77 |
|
|
æ -2.2531325651061 |
ö |
Находим вектор невязки: |
X_Jacobi |
= ç-2.40844472545144 |
÷ |
|
k |
ç |
÷ |
V_Jacobi := A×X_Jacobik - B |
|
è 6.59387169371675 |
ø |
|
æ |
− 9 |
ö |
ç |
7.64759988669539 ´ 10 |
÷ |
− 7 |
||
V_Jacobi = ç |
3.34393483569784 ´ 10 |
÷ |
ç |
− 7 |
÷ |
è |
-6.97594146004121 ´ 10 |
ø |
Norm1(V_Jacobi) = 1.0396352294606 ´ 10− 6
Norm2(V_Jacobi) = 7.73637563835406 ´ 10− 7
NormInf(V_Jacobi) = 6.97594146004121 ´ 10− 7
X_Zeidel := ZeidelPlus(A , B , 0 , e)
k2 := ZeidelPlus(A , B , 1 , e)
X - двумерный массив, где каждый из его элементов Xk- вектор-столбец, содержащий значения переменных системы, вычисленных на k-ом шаге.
Количество итераций k2 = 39
Находим значение вектора невязки: |
|
|
|
|
|
æ |
1.83383974672324 |
− 8 |
ö |
V_Zeidel := A×X_Zeidelk2 - B |
ç |
´ 10 |
÷ |
|
V_Zeidel = ç |
|
− 7 |
÷ |
|
|
ç |
3.54887042419705 |
´ 10 |
÷ |
|
è |
0 |
|
ø |
Norm1(V_Zeidel) = 3.73225439886937 ´ 10− 7
Norm2(V_Zeidel) = 3.55360534808062 ´ 10− 7
NormInf(V_Zeidel) = 3.54887042419705 ´ 10− 7
Функция нахождения максимальной из норм:
max_P(V) := s1 ¬ Norm1(V) s2 ¬ Norm2(V) s3 ¬ NormInf(V)
max(s)
Распечатаем промежуточные значения вектора-столбца переменных X, получаемых на различных промежуточных стадиях итерационных процессов Якоби и Зейделя:
Метод Якоби (77 итераций) |
Метод Зейделя (39 итераций) |
|
||||||
|
|
|
|
k := 2 |
|
|
|
|
|
æ -1.64738562091503 |
ö |
|
|
æ -1.64738562091503 |
ö |
||
X_Jacobi |
= ç |
-0.709313725490196 |
÷ |
X_Zeidel |
k |
= ç |
-0.936322568242984 |
÷ |
k |
ç |
|
÷ |
|
ç |
|
÷ |
|
|
è |
2.88848039215686 |
ø |
|
|
è |
5.18960367807254 |
ø |
|
|
|
|
|
k := 4 |
|
|
|
|
|
|
|
æ−1.92702081731385 ö |
|
|
|
æ−2.1671372712181 ö |
|
|||
X_Jacobi |
= |
ç |
−1.27548106497501 |
÷ |
X_Zeidel |
k |
= |
ç |
−1.8392645713063 ÷ |
|
k |
|
ç |
|
÷ |
|
|
ç |
÷ |
|
|
|
|
è 4.38042691913367 |
ø |
|
|
|
è |
6.05833713011345 ø |
|
|
|
|
|
|
|
k := 6 |
|
|
|
|
|
|
|
æ−2.07454154710952 ö |
|
|
|
æ−2.22039116976259 ö |
||||
X_Jacobi |
= |
ç |
−1.66445510813699 |
÷ |
X_Zeidel |
k |
= |
ç |
−2.19155667544109 |
÷ |
k |
|
ç |
|
÷ |
|
|
ç |
|
÷ |
|
|
|
è 5.26331138283761 |
ø |
|
|
|
è |
6.38980578889981 |
ø |
|
|
|
|
|
|
k := 8 |
|
|
|
|
|
|
|
æ−2.15452349870072 ö |
|
|
|
æ−2.24065643825257 ö |
||||
X_Jacobi |
= |
ç |
−1.925306645331 |
÷ |
X_Zeidel |
k |
= |
ç |
−2.32579939587395 |
÷ |
k |
|
ç |
|
÷ |
|
|
ç |
|
÷ |
|
|
|
è 5.78989808155743 |
ø |
|
|
|
è |
6.51611225806165 |
ø |
|
|
|
|
|
|
k := 10 |
|
|
|
|
|
|
|
æ−2.19821846201509 ö |
|
|
|
æ−2.24837852928473 ö |
||||
X_Jacobi |
= |
ç |
−2.09734938654213 |
÷ |
X_Zeidel |
k |
= |
ç |
−2.37695268794057 |
÷ |
k |
|
ç |
|
÷ |
|
|
ç |
|
÷ |
|
|
|
è 6.10587380863959 |
ø |
|
|
|
è |
6.56424143520877 |
ø |
|
|
|
|
|
|
k := 12 |
|
|
|
|
|
|
|
æ−2.22228019915953 ö |
|
|
|
æ−2.25132103800546 ö |
||||
X_Jacobi |
= |
ç |
−2.20943753174664 |
÷ |
X_Zeidel |
k |
= |
ç |
−2.3964446891903 |
÷ |
k |
|
ç |
|
÷ |
|
|
ç |
|
÷ |
|
|
|
è 6.29649332501364 |
ø |
|
|
|
è |
6.58258109514456 |
ø |
|
|
|
|
|
|
k := 14 |
|
|
|
|
|
|
|
æ−2.23564240331273 ö |
|
|
|
æ −2.2524422832146 |
ö |
|||
X_Jacobi |
= |
ç |
−2.28179238556335 |
÷ |
X_Zeidel |
k |
= |
ç |
−2.40387213124474 |
÷ |
k |
|
ç |
|
÷ |
|
|
ç |
|
÷ |
|
|
|
è 6.41203423288301 |
ø |
|
|
|
è |
6.58956943665582 |
ø |
|
|
|
|
|
|
k := 16 |
|
|
|
|
|
|
|
æ−2.24312809809888 ö |
|
|
|
æ−2.25286953456314 ö |
||||
X_Jacobi |
= |
ç |
−2.32816770306463 |
÷ |
X_Zeidel |
k |
= |
ç |
−2.40670236374705 |
÷ |
k |
|
ç |
|
÷ |
|
|
ç |
|
÷ |
|
|
|
è 6.48235754594489 |
ø |
|
|
|
è |
6.59223234955874 |
ø |
|
|
|
|
|
k := 18 |
|
|
|
|
|
|
æ-2.24735939915543 |
ö |
|
|
|
æ-2.25303233901782 ö |
|
X_Jacobi |
= |
ç-2.35772708809334 |
÷ |
X_Zeidel |
k |
= |
ç-2.40778082601689 |
÷ |
k |
|
ç |
÷ |
|
|
ç |
÷ |
|
|
|
è 6.52531299160262 |
ø |
|
|
|
è 6.59324705456666 |
ø |
Построим график, отражающий динамику сходимости методов Якоби и Зейделя:
get_Conv_graph(A , X , B) := for k Î 2 .. length(X)
Convk−1 ¬ max_P(A×Xk - B)
Conv
Jacobi_conv := get_Conv_graph(A , X_Jacobi, B)
Zeidel_conv := get_Conv_graph(A , X_Zeidel , B)
Jacobi_conv |
Zeidel_conv |
Произведем сравнение использованных методов по векторам невязки:
Norm1(V_Gauss) = 2.1316282072803 ´ 10− 14
Norm2(V_Gauss) = 1.50728876033642 |
´ 10− 14 |
|
|
NormInf(V_Gauss) = 1.4210854715202 |
´ 10− |
14 |
|
Norm1(V_Jacobi) = 1.0396352294606 ´ 10− 6 |
Norm1(V_Zeidel) = 3.73225439886937 ´ 10− 7 |
||
Norm2(V_Jacobi) = 7.73637563835406 |
´ 10− |
7 |
Norm2(V_Zeidel) = 3.55360534808062 ´ 10− 7 |
NormInf(V_Jacobi) = 6.97594146004121 ´ 10− 7 |
NormInf(V_Zeidel) = 3.54887042419705 ´ 10− 7 |
||
Вывод: Как видно из вышепредставленных промежуточных значений итерационного процесса, метод Зейделя сходится быстрее (примерно в два раза), хотя оба метода имеют геометрическую сходимость. Наименьшая погрешность была достигнута при использовании метода Гаусса с целочисленной арифметикой. В
ходе вычислительного эксперимента метод Зейделя показывал тем лучшие результаты вектора невязки, чем ниже в векторе-столбце находилась переменная. Это связано с тем, что каждая из последующих переменных вычисляется уже на основе только что вычисленных предыдущих, что обеспечивает более глубокое автокорректирование результата. И когда норма, согласованная с нормой матрицы P, будет отражать факт достижения заданной точности, точность вычисления последней переменной в векторе-столбце X будет тем значительнее превосходить заданную, чем больше уравнений было в исходной системе. Целочисленный метод Гаусса превосходит использованные итерационные методы по точности в два раза (при условии представления конечного результата целочисленного метода Гаусса в вещественном виде).
При представлении же результата целочисленного метода Гаусса в натуральных дробях погрешность вообще отсутсвует как таковая.
Задача 2.7. Дана система уравнений AХ=В порядка n. Исследовать зависимость погрешности решения Х от
погрешностей элементов правой части системы В. ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:
2.Задать матрицу системы A и вектор правой части В. Используя программу для решения систем методом Гаусса с выбором ведущего элемента найти решение системы.
3.Принять решение Х, полученное в п. 1, за точное Х*. Поочередно изменяя значение m-го элемента
вектор - столбца свободных членов В на 1% от исходного значения в сторону увеличения по модулю, находить решение Хвозм(m) возмущенной системы методом Гаусса с выбором ведущего элемента.
4.Составить вектор отклонений d = (d1 ,d2 ,...,dn ) по формуле
dm = |
X * − X возм(m) |
C |
|||
|
|
|
|||
|
|
|
|
|
|
|
|
|
X * |
C |
|
|
|
|
|
||
На основе вычисленного вектора d построить гистограмму. По гистограмме определить компоненту bm вектора В, которая оказывает наибольшее влияние на погрешность решения. Сделать выводы
Инициализируем исходную систему:
n := 100 + N |
n = 114 |
||
fill_A := |
|
for |
i Î 1 .. n |
|
|||
|
|
Ai, 1 ¬ 1 |
|
|
|
for |
k Î 1 .. n - 1 |
|
|
|
Ak+1, k+1 ¬ 1 + (k + N) ×0.01 |
|
|
|
|
|
|
|
Ak+2, k+1 ¬ 1 if k < n - 1 |
|
|
A1, n ¬ 1 |
|
|
|
A |
|
|
|
|
1 |
|
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
|
|
1 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
2 |
|
1 |
1.15 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
3 |
|
1 |
1 |
1.16 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
4 |
|
1 |
0 |
1 |
1.17 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
5 |
|
1 |
0 |
0 |
1 |
1.18 |
0 |
0 |
0 |
0 |
0 |
|
|
6 |
|
1 |
0 |
0 |
0 |
1 |
1.19 |
0 |
0 |
0 |
0 |
A := fill_A |
A = |
7 |
|
1 |
0 |
0 |
0 |
0 |
1 |
1.2 |
0 |
0 |
0 |
8 |
|
1 |
0 |
0 |
0 |
0 |
0 |
1 |
1.21 |
0 |
0 |
||
|
|
9 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
1.22 |
0 |
|
|
10 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
1.23 |
|
|
11 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
|
|
12 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
13 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
14 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
15 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
16 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
fill_B := |
|
for |
k 1 .. n − 1 |
|
|
|
|||||
|
|
|
1 |
||||||||
|
|
|
|
é |
(π + N) ù |
|
1 |
0.149798873510009 |
|||
|
|
Bk ← sinê |
|
|
ú |
|
|
|
|||
|
|
|
|
|
2 |
0.151114388328177 |
|||||
|
|
|
|
|
ën − k + 1û |
|
|||||
|
|
|
← sin |
é |
(π + N) |
ù |
|
3 |
0.152453119685847 |
||
|
|
Bn |
|
4 |
0.153815685040195 |
||||||
|
|
ë |
|
1 û |
|
||||||
|
|
B |
|
|
|
|
5 |
0.155202723840962 |
|||
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
6 |
0.156614898513671 |
|
|
|
|
|
|
|
|
|
|
|
||
B := fill_B |
|
|
|
|
|
|
|
B = |
7 |
0.158052895495781 |
|
|
|
|
|
|
|
|
8 |
0.159517426329128 |
|||
|
|
|
|
|
|
|
|
|
|
9 |
0.161009228812226 |
|
|
|
|
|
|
|
|
|
|
10 |
0.162529068216261 |
|
|
|
|
|
|
|
|
|
|
11 |
0.164077738568911 |
|
|
|
|
|
|
|
|
|
|
12 |
0.165656064010406 |
|
|
|
|
|
|
|
|
|
|
13 |
0.167264900226587 |
|
|
|
|
|
|
|
|
|
|
14 |
0.168905135964086 |
|
|
|
|
|
|
|
|
|
|
15 |
0.170577694633126 |
|
|
|
|
|
|
|
|
|
|
16 |
0.172283536003865 |
|
|
|
|
|
|
|
|
|
|
|
|
Сначала найдём решение невозмущенной системы: |
|
|
|
|
||||||||
|
1 |
|
|
|||||||||
X_true := |
GaussReverse(GaussForwardPlus(A , B)) |
1 |
1.08047963668287 |
|
||||||||
2 |
-0.808143694221472 |
|
||||||||||
|
|
|
|
|
|
|
|
|
3 |
-0.103347261013406 |
|
|
|
|
|
|
|
|
|
|
|
4 |
-0.703689479170316 |
|
|
|
|
|
|
|
|
|
|
|
5 |
-0.187785960738637 |
|
|
|
|
|
|
|
|
|
|
|
6 |
-0.618553594479463 |
|
|
|
|
|
|
|
|
|
|
X_true = |
7 |
-0.253227622256355 |
|
|
|
|
|
|
|
|
|
|
8 |
-0.551846767022634 |
|
||
|
|
|
|
|
|
|
|
|
9 |
-0.301330853154107 |
|
|
|
|
|
|
|
|
|
|
|
10 |
-0.501316841717482 |
|
|
|
|
|
|
|
|
|
|
|
11 |
-0.334746013222965 |
|
|
|
|
|
|
|
|
|
|
|
12 |
-0.464062047559599 |
|
|
|
|
|
|
|
|
|
|
|
13 |
-0.356470388013242 |
|
|
|
|
|
|
|
|
|
|
|
14 |
-0.43708985252405 |
|
|
|
|
|
|
|
|
|
|
|
15 |
-0.369384444941949 |
|
|
|
|
|
|
|
|
|
|
|
16 |
-0.417683454059734 |
|
|
get_D := |
|
C ← B |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
for i 1 .. n |
|
|
|
|
|
|
|
|||
|
|
|
Ci |
← Ci |
+ |
Bi |
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
100 |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
Di |
← |
NormInf((X_true − GaussReverse(GaussForwardPlus(A , C)))) |
D := get_D |
||||||
|
|
|||||||||||
|
|
|
|
|
|
|
|
NormInf(X_true) |
|
|
|
|
|
|
D |
|
|
|
|
|
|
|
|
|
|
Построим гистограмму погрешностей результата в зависимости от номера элемента вектора-столбца свободных членов.
D
Вывод: Как видно из построенной гистограммы, наибольшее влияние на отклонение
от результата при одинаковом процентном возмущении получается при внесении погрешности в наибольший по модулю элемент вектора-столбца свободных членов. Это объясняется тем, что наибольший по модулю элемент
при равных процентных возмущениях вносит наибольшую абсолютную погрешность в исходные данные.