Z-ট্রান্সফর্ম দিয়ে ডিফারেন্স ইকুয়েশন সমাধান

import cmath

# ডিফারেন্স ইকুয়েশন: y[n] - 1.5 y[n-1] + 0.5 y[n-2] = 0   (zero-input response)
# প্রাথমিক শর্ত: y[0] = 2, y[1] = 0.5

print("=== ধাপ ১: characteristic equation ও পোল (L18-এর পদ্ধতি পুনর্ব্যবহার) ===")
a, b, c = 1, -1.5, 0.5
disc = b*b - 4*a*c
print(f"z^2 + ({b})z + ({c}) = 0,  ডিসক্রিমিন্যান্ট = {disc}")
sqrt_disc = cmath.sqrt(disc)
z1 = (-b + sqrt_disc) / (2*a)
z2 = (-b - sqrt_disc) / (2*a)
print(f"পোল ১: z1 = {z1}")
print(f"পোল ২: z2 = {z2}")

print()
print("=== ধাপ ২: সাধারণ সমাধান y[n] = C1*z1^n + C2*z2^n, প্রাথমিক শর্ত থেকে C1, C2 নির্ণয় ===")
y0, y1 = 2, 0.5
# C1 + C2 = y0
# C1*z1 + C2*z2 = y1
C1 = (y1 - y0*z2) / (z1 - z2)
C2 = (y0*z1 - y1) / (z1 - z2)
print(f"C1 = {C1}")
print(f"C2 = {C2}")

def y_closed(n):
    return (C1 * z1**n + C2 * z2**n).real

print()
print("=== ধাপ ৩: ক্লোজড-ফর্ম বনাম সরাসরি সিমুলেশন -- ক্রস-চেক ===")
N = 10
y_sim = [0.0]*N
y_sim[0] = y0
y_sim[1] = y1
for n in range(2, N):
    y_sim[n] = 1.5*y_sim[n-1] - 0.5*y_sim[n-2]

print(f"{'n':>3} | {'closed-form':>14} | {'সরাসরি সিমুলেশন':>16} | {'পার্থক্য':>10}")
for n in range(N):
    cf = y_closed(n)
    diff = abs(cf - y_sim[n])
    print(f"{n:>3} | {cf:>14.8f} | {y_sim[n]:>16.8f} | {diff:>10.2e}")

print()
print(f"পোল z1 = {z1.real} (|z1|={abs(z1):.4f}, ইউনিট সার্কেলের উপর) -> এই মোড কমে না, ধ্রুবক থাকে")
print(f"পোল z2 = {z2.real} (|z2|={abs(z2):.4f} < 1) -> এই মোড ক্ষয়প্রাপ্ত হয়")