Метод простой итерации и метод Зейделя

@VictorVAlduin · Регистрация: 15.12.2020

Студворк — интернет-сервис помощи студентам

Нужна помощь плз!
Есть ли у кого-то ресурсы пример реализации какого-либо из методов с помощью библиотеки numpy
Погулил.Все что нашел реализации данных методов вручную "пальчиками" и код ужасно огромен.
Вот например:
# Метод простых итераций и метод Зейделя
# Вариант

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
Nc = 12
U = Nc + 10
 
f1 = [ 0., (-1./U), (-1./U), (1./U) ]
f2 = [ (-1./U), 0, (-3./U), 1. ]
f3 = [ (-1./5), (-3./5), 0, 0 ]
###########
# Simple x,y,z
def x(n):
    global GT
    if len(GT[0])>n:
        return GT[0][n]
    else:
        return f1[3] if n<=0 else f1[0]*x(n-1)+f1[1]*y(n-1)+f1[2]*z(n-1)+f1[3]
def y(n):
    global GT
    if len(GT[1])>n:
        return GT[1][n]
    else:
        return f2[3] if n<=0 else f2[0]*x(n-1)+f2[1]*y(n-1)+f2[2]*z(n-1)+f2[3]
def z(n):
    global GT
    if len(GT[2])>n:
        return GT[2][n]
    else:
        return f3[3] if n<=0 else f3[0]*x(n-1)+f3[1]*y(n-1)+f3[2]*z(n-1)+f3[3]
# Seidel x,y,z
def Sx(n):
    global ST
    if len(ST[0])>n:
        return ST[0][n]
    else:
        return 0 if n<=0 else f1[0]*Sx(n-1)+f1[1]*Sy(n-1)+f1[2]*Sz(n-1)+f1[3]
def Sy(n):
    global ST
    if len(ST[1])>n:
        return ST[1][n]
    else:
        return 0 if n<=0 else f2[0]*Sx(n)+f2[1]*Sy(n-1)+f2[2]*Sz(n-1)+f2[3]
def Sz(n):
    global ST
    if len(ST[2])>n:
        return ST[2][n]
    else:
        return 0 if n<=0 else f3[0]*Sx(n)+f3[1]*Sy(n)+f3[2]*Sz(n-1)+f3[3]
###########
lBl1 = max( abs(f1[0])+abs(f1[1])+abs(f1[2]) , abs(f2[0])+abs(f2[1])+abs(f2[2]) , abs(f3[0])+abs(f3[1])+abs(f3[2]) )
lBl2 = max( abs(f1[0])+abs(f2[0])+abs(f3[0]) , abs(f1[1])+abs(f2[1])+abs(f3[1]) , abs(f1[2])+abs(f2[2])+abs(f3[2]) )
lBl3 = ( abs(f1[0])**2+abs(f2[0])**2+abs(f3[0])**2+abs(f1[1])**2+abs(f2[1])**2+abs(f3[1])+abs(f1[2])**2+abs(f2[2])**2+abs(f3[2])**2 )**(0.5)
lBl = min( lBl1 , lBl2, lBl3)
 
if lBl>=1:
    print "Итерационный процесс не сходится к точному решению"
    exit()
    
if lBl==lBl1:
    ldl = max( abs(f1[3]), abs(f2[3]), abs(f3[3]) )
elif lBl==lBl2:
    ldl = abs(f1[3]) + abs(f2[3]) + abs(f3[3])
else:
    ldl = ( abs(f1[3])**2 + abs(f2[3])**2 + abs(f3[3])**2 )**(0.5)
 
def E(n):
    global lBl, ldl
    return (lBl**n * ldl) / ( 1 - lBl )
 
delta = lambda n: max( abs(x(n)-x(n-1)), abs(y(n)-y(n-1)), abs(z(n)-z(n-1)) )
Sdelta = lambda n: max( abs(Sx(n)-Sx(n-1)), abs(Sy(n)-Sy(n-1)), abs(Sz(n)-Sz(n-1)) )
GT = [ [], [], [] ]
ST = [ [], [], [] ]
 
GT[0].append(x(0))
GT[1].append(y(0))
GT[2].append(z(0))
 
print "Simple iteration method. Variant: %i" % Nc
print "%2i: %2.16f %2.16f %2.16f" % (0, x(0), y(0), z(0))
i = 1
while ( E(i-1)>10**(-5) ):
    GT[0].append(x(i))
    GT[1].append(y(i))
    GT[2].append(z(i))
    print "%2i: %2.16f %2.16f %2.16f %2.16f %2.16f" % (i, x(i), y(i), z(i), E(i), delta(i))
    i += 1
 
print "\n\r\n\rSeidel method. Variant: %i" % Nc
print "%2i: %2.16f %2.16f %2.16f" % (0, Sx(0), Sy(0), Sz(0))
i = 1
while ( Sdelta(i)>0 ):
    ST[0].append(Sx(i))
    ST[1].append(Sy(i))
    ST[2].append(Sz(i))
    print "%2i: %2.16f %2.16f %2.16f %2.16f %2.16f" % (i, Sx(i), Sy(i), Sz(i), E(i), Sdelta(i))
    i += 1

Поэтому спрашиваю есть ли реализацию через библиотеки что упростило б жизнь.Или и вправду нужно ручками будет все набирать
Буду благодарен за представленные ресурсы если есть!
Для каждого метода нужно вывести количество
итераций, необходимых для разрешения СЛАР.

@Aael · 27.09.2022, 23:48

Сообщение от VictorVAlduin

Поэтому спрашиваю есть ли реализацию через библиотеки что упростило б жизнь.Или и вправду нужно ручками будет все набирать!

Так в серьезных либах не используют классические методы из учебников. Но если они тебе так нужны, то почему бы и не набрать пальчиками? Тем более весь метод Зейделя к микроформуле сводится (с numpy тоже, матричной арифметикой):

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Ax = b. Метод Зейделя
A = [[-0.76,-0.04,0.21,-0.18],[0.45,-1.23,0.66,0.0],[0.26,0.34,-1.11,0.0],[0.05,-0.26,0.34,-1.12]]
b = [-1.24,0.88,-0.63,1.17]
 
from math import dist
x, y, n = b, [0] * len(b), 0
while dist(x, y) > 1e-6:
    x, y, n = [0] * len(x), x, n+1
    for i in range(len(y)):
        x[i] = 1 / A[i][i] * (b[i] -
            sum(x[j] * A[i][j] for j in range(i)) - 
            sum(y[j] * A[i][j] for j in range(i+1, len(y))))
 
# x = (2.12, 0.76, 1.30, -0.73) итераций: 12
print(f'x = ({", ".join(f"{xi:.2f}" for xi in x)}) итераций: {n}')

@VictorVAlduin · 28.09.2022, 00:09 **[ТС]**

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def seidel(A, b, eps):
    n = len(A)
    x = np.zeros(n)  # zero vector
 
    converge = False
    while not converge:
        x_new = np.copy(x)
        for i in range(n):
            s1 = sum(A[i][j] * x_new[j] for j in range(i))
            s2 = sum(A[i][j] * x[j] for j in range(i + 1, n))
            x_new[i] = (b[i] - s1 - s2) / A[i][i]
 
        converge = np.sqrt(sum((x_new[i] - x[i]) ** 2 for i in range(n))) <= eps
        x = x_new
 
    return x
 
A = [[-0.76, -0.04, 0.21, -0.18], [0.45, -1.23, 0.66, 0.0], [0.26, 0.34, -1.11, 0.0], [0.05, -0.26, 0.34, -1.12]]
b = [-1.24, 0.88, -0.63, 1.17]
seidel(A, b, 1e-4))

Добавлено через 3 минуты
Хм у меня вышло 4 итерации.Это нормально ?

@Aael · 28.09.2022, 01:10

VictorVAlduin, у вас другое начальное приближение (нули, у меня - b) и искомая точность разная. Матричная формула для нумпая (норма с abs вместо sqrt, можно поменять):

Python
1
2
3
4
5
6
7
8
9
10
11
import numpy as np
     
A = np.array([[-0.76, -0.04, 0.21, -0.18], [0.45, -1.23, 0.66, 0.0], [0.26, 0.34, -1.11, 0.0], [0.05, -0.26, 0.34, -1.12]])
b = np.array([-1.24, 0.88, -0.63, 1.17])
 
x, y, n = b.copy(), np.zeros(b.size), 0
U, L = -np.triu(A, 1), np.linalg.inv(np.tril(A))
while np.abs(x-y).max() > 1e-6:
    x, y, n = L.dot(U.dot(x) + b), x, n+1
 
print(f'{x=} итераций={n}')

@VictorVAlduin · 28.09.2022, 02:39 **[ТС]**

Я хотел вывести числа с точностью до 4 знака .Но выводится почему то больше.
А также тут я сделал так чтобы типо я определяю сколько будет количество итерации.
Стоит ли так оставлять или лучше модернизировать и сделать чтобы оно само считало итерации?Если да подскажите опять же как тут лучше поступить?
Это метод простой итерации (метод Якоби)

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def jacobi(A,b,N=25,x=None):
    """Solves the equation Ax=b via the Jacobi iterative method."""
    # Create an initial guess if needed
    if x is None:
        x = np.zeros(len(A[0]))
 
    # Create a vector of the diagonal elements of A
    # and subtract them from A
    D = np.diag(A)
    R = A - np.diagflat(D)
 
    # Iterate for N times
    for i in range(N):
        x = (b - np.dot(R,x)) / D
    return x
 
A = np.array([[-0.76, -0.04, 0.21, -0.18], [0.45, -1.23, 0.66, 0.0],
              [0.26, 0.34, -1.11, 0.0], [0.05, -0.26, 0.34, -1.12]])
b = np.array([-1.24, 0.88, -0.63, 1.17])
guess = np.array([1e-4,1e-4,1e-4,1e-4])
 
sol = jacobi(A,b,N=25,x=guess)
 
print("A:")
pprint(A)
 
print("b:")
pprint(b)
 
print("x:")
pprint(sol)

Добавлено через 1 час 1 минуту
Проблема решена.Все получилось

@u235 · 28.09.2022, 08:14

Python
1
2
3
4
5
6
7
8
9
10
import numpy as np
a=np.array([[-.76, -.04,.21,-.18], [.45,-1.23,.66, 0],[.26,.34,-1.11, 0],[.05,-.26,.34,-1.12]])
b=np.array([-1.24, .88, -.63, 1.17])
alpha=-a/np.tile(np.diag(a), (4,1))
alpha[np.diag_indices_from(alpha)]=0
beta=b/np.diag(a)
x=beta
for i in range(10):
    x=alpha@x+beta
    print(x)

Но где-то у меня ошибка, в том смысле что решение не сходится к правильному.. Ищите в чем причина.

@VictorVAlduin · 28.09.2022, 14:32 **[ТС]**

Что я делаю не так ?

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# Defining equations to be solved
# in diagonally dominant form
f1 = lambda x, y, z, d: (1.24 - 0.04*y + 0.21 * z - 0.18 * d) / 0.76
f2 = lambda x, y, z, d: (-0.88 + 0.45 * x + 0.66 * z) / 1.23
f3 = lambda x, y, z, d: (0.63 + 0.26 * x + 0.34 * y) / 1.11
f4 = lambda x, y, z, d: (-1.17 + 0.05 * x - 0.26 * y + 0.34 * z) / 1.12
 
# Initial setup
x0 = 0
y0 = 0
z0 = 0
d0 = 0
count = 1
 
# Reading tolerable error
e = float(input('Enter tolerable error: '))
 
# Implementation of Jacobi Iteration
print('\nCount\tx\ty\tz\td\n')
 
condition = True
 
while condition:
    x1 = f1(x0, y0, z0, d0)
    y1 = f2(x0, y0, z0, d0)
    z1 = f3(x0, y0, z0, d0)
    d1 = f4(x0, y0, z0, d0)
    print('%d\t%0.4f\t%0.4f\t%0.4f\n' % (count, x1, y1, z1, d1))
    e1 = abs(x0 - x1)
    e2 = abs(y0 - y1)
    e3 = abs(z0 - z1)
    e4 = abs(d0 - d1)
 
    count += 1
    x0 = x1
    y0 = y1
    z0 = z1
    d0 = d1
 
    condition = e1 > e and e2 > e and e3 > e and e4 > e
 
print('\nSolution: x=%0.4f, y=%0.4f and z = %0.4f\n' % (x1, y1, z1, d1))

Подскажите как правильно сделать?Спасибо заранее
Хочу реализовать метод за формулами:

@Aael · 28.09.2022, 14:49

VictorVAlduin, нехватает еще одного %0.4f в 28 и 42 строках для компоненты d

@VictorVAlduin · 28.09.2022, 17:29 **[ТС]**

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def simple_iter(f1, f2, f3, f4):
 
    x0 = 0
    y0 = 0
    z0 = 0
    d0 = 0
    count = 1
    n = 0
 
    e = float(input('Точность: '))
 
 
    print('\nCount\tx\ty\tz\td\n')
 
    condition = True
 
    while condition:
        x1 = f1(x0, y0, z0, d0)
        y1 = f2(x0, y0, z0, d0)
        z1 = f3(x0, y0, z0, d0)
        d1 = f4(x0, y0, z0, d0)
        print('%d\t%0.4f\t%0.4f\t%0.4f \t%0.4f\n' % (count, x1, y1, z1, d1))
        e1 = abs(x0 - x1)
        e2 = abs(y0 - y1)
        e3 = abs(z0 - z1)
        e4 = abs(d0 - d1)
        n += 1
        count += 1
        x0 = x1
        y0 = y1
        z0 = z1
        d0 = d1
 
        condition = e1 > e and e2 > e and e3 > e and e4 > e
 
    print('\nРешения: x=%0.4f, y=%0.4f and z = %0.4f d = %0.4f\n' % (x1, y1, z1, d1), f'\nИтерации: {n=}')
 
 
f1 = lambda x, y, z, d: (1.24 - 0.04*y + 0.21 * z - 0.18 * d) / 0.76
f2 = lambda x, y, z, d: (-0.88 + 0.45 * x + 0.66 * z) / 1.23
f3 = lambda x, y, z, d: (0.63 + 0.26 * x + 0.34 * y) / 1.11
f4 = lambda x, y, z, d: (-1.17 + 0.05 * x - 0.26 * y + 0.34 * z) / 1.12

Тут где-то ошибка в плане итерации выходит 8 а должно 15.
Что не так делаю?

Новые блоги и статьи Все статьи Все блоги /
Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Programma_Boinc 28.12.2025 Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Налог на собак: https:/ / ********/ gallery/ V06K53e Финансовый отчет в Excel: https:/ / ********/ gallery/ bKBkQFf Пост отсюда. . .	Кто-нибудь знает, где можно бесплатно получить настольный компьютер или ноутбук? США. Programma_Boinc 26.12.2025 Нашел на реддите интересную статью под названием Anyone know where to get a free Desktop or Laptop? Ниже её машинный перевод. После долгих разбирательств я наконец-то вернула себе. . .	Thinkpad X220 Tablet — это лучший бюджетный ноутбук для учёбы, точка. Programma_Boinc 23.12.2025 Рецензия / Мнение/ Перевод Нашел на реддите интересную статью под названием The Thinkpad X220 Tablet is the best budget school laptop period . Ниже её машинный перевод. Thinkpad X220 Tablet —. . .	PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Как объединить две одинаковые БД Access с разными данными VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.
Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .

@Aael 964 / 485 / 241 Регистрация: 02.06.2016 Сообщений: 760
	28.09.2022, 14:49
	VictorVAlduin, нехватает еще одного `%0.4f` в 28 и 42 строках для компоненты d 1