ইটারেটিভ মেথড — জ্যাকোবি ও গস-সেইডেল

def jacobi_step(A, b, x_old):
    n = len(A)
    x_new = [0.0] * n
    for i in range(n):
        s = sum(A[i][j] * x_old[j] for j in range(n) if j != i)
        x_new[i] = (b[i] - s) / A[i][i]
    return x_new

def gauss_seidel_step(A, b, x):
    n = len(A)
    x_new = x[:]
    for i in range(n):
        s = sum(A[i][j] * x_new[j] for j in range(n) if j != i)
        x_new[i] = (b[i] - s) / A[i][i]
    return x_new

def run(A, b, method, n_iter, cap=1e12):
    n = len(A)
    x = [0.0] * n
    history = [x[:]]
    for k in range(n_iter):
        x = jacobi_step(A, b, x) if method == "jacobi" else gauss_seidel_step(A, b, x)
        history.append(x[:])
        if any(abs(v) > cap for v in x):
            break   # ব্লো-আপ থামিয়ে দেওয়া
    return history

# --- ডায়াগোনালি-ডমিন্যান্ট সিস্টেম (কনভার্জ করবে) ---
A_conv = [
    [10.0, -1.0,  2.0],
    [ 1.0, 11.0, -1.0],
    [ 2.0, -1.0, 10.0],
]
b_conv = [6.0, 25.0, -11.0]

hist_j = run(A_conv, b_conv, "jacobi", 10)
hist_g = run(A_conv, b_conv, "gauss_seidel", 10)

print("Jacobi ইটারেশন টেবিল:")
print(f"{'k':>2} | {'x1':>10} {'x2':>10} {'x3':>10}")
for k, xk in enumerate(hist_j):
    print(f"{k:>2} | {xk[0]:>10.6f} {xk[1]:>10.6f} {xk[2]:>10.6f}")

print()
print("Gauss-Seidel ইটারেশন টেবিল:")
print(f"{'k':>2} | {'x1':>10} {'x2':>10} {'x3':>10}")
for k, xk in enumerate(hist_g):
    print(f"{k:>2} | {xk[0]:>10.6f} {xk[1]:>10.6f} {xk[2]:>10.6f}")

# L10-এর গসিয়ান এলিমিনেশন দিয়ে রেফারেন্স সমাধান
def gaussian_solve(A, b):
    n = len(A)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for k in range(n - 1):
        for i in range(k + 1, n):
            f = M[i][k] / M[k][k]
            for j in range(k, n + 1):
                M[i][j] -= f * M[k][j]
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        s = M[i][n]
        for j in range(i + 1, n):
            s -= M[i][j] * x[j]
        x[i] = s / M[i][i]
    return x

x_ref = gaussian_solve(A_conv, b_conv)
print()
print("রেফারেন্স সমাধান (Gaussian elimination):", [round(v, 6) for v in x_ref])
err_j = max(abs(hist_j[-1][i] - x_ref[i]) for i in range(3))
err_g = max(abs(hist_g[-1][i] - x_ref[i]) for i in range(3))
print(f"১০ ইটারেশন পর Jacobi সর্বোচ্চ ভুল: {err_j:.2e}")
print(f"১০ ইটারেশন পর Gauss-Seidel সর্বোচ্চ ভুল: {err_g:.2e}")

# --- ডায়াগোনালি-ডমিন্যান্ট নয় এমন সিস্টেম (ডাইভার্জ করবে) ---
print()
print("=== ডায়াগোনালি-ডমিন্যান্ট নয় এমন সিস্টেম ===")
A_div = [[1.0, 2.0], [3.0, 1.0]]
b_div = [5.0, 10.0]

hist_jd = run(A_div, b_div, "jacobi", 10)
print("Jacobi ইটারেশন টেবিল (ডাইভার্জিং):")
print(f"{'k':>2} | {'x1':>16} {'x2':>16}")
for k, xk in enumerate(hist_jd):
    print(f"{k:>2} | {xk[0]:>16.4f} {xk[1]:>16.4f}")

Jacobi ইটারেশন টেবিল:
 k |         x1         x2         x3
 0 |   0.000000   0.000000   0.000000
 1 |   0.600000   2.272727  -1.100000
 2 |   1.047273   2.118182  -0.992727
 3 |   1.010364   2.087273  -1.097636
 4 |   1.028255   2.081091  -1.093345
 5 |   1.026778   2.079855  -1.097542
 6 |   1.027494   2.079607  -1.097370
 7 |   1.027435   2.079558  -1.097538
 8 |   1.027463   2.079548  -1.097531
 9 |   1.027461   2.079546  -1.097538
10 |   1.027462   2.079546  -1.097538

Gauss-Seidel ইটারেশন টেবিল:
 k |         x1         x2         x3
 0 |   0.000000   0.000000   0.000000
 1 |   0.600000   2.218182  -0.998182
 2 |   1.021455   2.089124  -1.095379
 3 |   1.027988   2.079694  -1.097628
 4 |   1.027495   2.079534  -1.097546
 5 |   1.027463   2.079545  -1.097538
 6 |   1.027462   2.079545  -1.097538
 7 |   1.027462   2.079545  -1.097538
 8 |   1.027462   2.079545  -1.097538
 9 |   1.027462   2.079545  -1.097538
10 |   1.027462   2.079545  -1.097538

রেফারেন্স সমাধান (Gaussian elimination): [1.027462, 2.079545, -1.097538]
১০ ইটারেশন পর Jacobi সর্বোচ্চ ভুল: 2.68e-07
১০ ইটারেশন পর Gauss-Seidel সর্বোচ্চ ভুল: 1.95e-13

=== ডায়াগোনালি-ডমিন্যান্ট নয় এমন সিস্টেম ===
Jacobi ইটারেশন টেবিল (ডাইভার্জিং):
 k |               x1               x2
 0 |           0.0000           0.0000
 1 |           5.0000          10.0000
 2 |         -15.0000          -5.0000
 3 |          15.0000          55.0000
 4 |        -105.0000         -35.0000
 5 |          75.0000         325.0000
 6 |        -645.0000        -215.0000
 7 |         435.0000        1945.0000
 8 |       -3885.0000       -1295.0000
 9 |        2595.0000       11665.0000
10 |      -23325.0000       -7775.0000