পিভোটিং স্ট্র্যাটেজি ও নিউমেরিক্যাল স্ট্যাবিলিটি

import math

def round_sig(x, sig=3):
    # সিমুলেট করে একটি সীমিত-নির্ভুলতার (limited-precision) মেশিন -- প্রতিটি arithmetic
    # অপারেশনের পর ফলাফল নির্দিষ্টসংখ্যক সিগনিফিক্যান্ট ডিজিটে রাউন্ড করা হয়
    if x == 0:
        return 0.0
    d = math.ceil(math.log10(abs(x)))
    power = sig - d
    factor = 10 ** power
    return round(x * factor) / factor

def solve_2x2(A, b, pivot, sig=3):
    M = [A[0][:] + [b[0]], A[1][:] + [b[1]]]

    if pivot and abs(M[1][0]) > abs(M[0][0]):
        M[0], M[1] = M[1], M[0]

    factor = round_sig(M[1][0] / M[0][0], sig)
    for j in range(3):
        M[1][j] = round_sig(M[1][j] - round_sig(factor * M[0][j], sig), sig)

    x2 = round_sig(M[1][2] / M[1][1], sig)
    x1 = round_sig(round_sig(M[0][2] - round_sig(M[0][1] * x2, sig), sig) / M[0][0], sig)
    return [x1, x2]

A = [
    [0.0001, 1.0],
    [1.0,    1.0],
]
b = [1.0, 2.0]

# পূর্ণ-নির্ভুলতায় (double precision) প্রকৃত সমাধান -- Cramer's rule দিয়ে রেফারেন্স
det = A[0][0] * A[1][1] - A[0][1] * A[1][0]
true_x1 = (b[0] * A[1][1] - A[0][1] * b[1]) / det
true_x2 = (A[0][0] * b[1] - b[0] * A[1][0]) / det
print(f"রেফারেন্স সমাধান (পূর্ণ ডাবল-প্রিসিশন): x1 = {true_x1:.8f}   x2 = {true_x2:.8f}")

sig = 3
x_naive = solve_2x2(A, b, pivot=False, sig=sig)
x_pivot = solve_2x2(A, b, pivot=True, sig=sig)

print()
print(f"{sig}-সিগনিফিক্যান্ট-ডিজিট মেশিনে, পিভোটিং ছাড়া (naive):")
print(f"  x1 = {x_naive[0]}   x2 = {x_naive[1]}")
print(f"  x1-এর প্রকৃত ভুল = {abs(x_naive[0] - true_x1):.6f}   x2-এর প্রকৃত ভুল = {abs(x_naive[1] - true_x2):.6f}")

print()
print(f"{sig}-সিগনিফিক্যান্ট-ডিজিট মেশিনে, partial pivoting সহ:")
print(f"  x1 = {x_pivot[0]}   x2 = {x_pivot[1]}")
print(f"  x1-এর প্রকৃত ভুল = {abs(x_pivot[0] - true_x1):.6f}   x2-এর প্রকৃত ভুল = {abs(x_pivot[1] - true_x2):.6f}")

রেফারেন্স সমাধান (পূর্ণ ডাবল-প্রিসিশন): x1 = 1.00010001   x2 = 0.99989999

3-সিগনিফিক্যান্ট-ডিজিট মেশিনে, পিভোটিং ছাড়া (naive):
  x1 = 0.0   x2 = 1.0
  x1-এর প্রকৃত ভুল = 1.000100   x2-এর প্রকৃত ভুল = 0.000100

3-সিগনিফিক্যান্ট-ডিজিট মেশিনে, partial pivoting সহ:
  x1 = 1.0   x2 = 1.0
  x1-এর প্রকৃত ভুল = 0.000100   x2-এর প্রকৃত ভুল = 0.000100