স্টিফ সমীকরণ ও ODE সলভারের স্ট্যাবিলিটি

১ · স্টিফ সমীকরণ কী

L27–L31 পর্যন্ত আমরা "নির্ভুলতা" নিয়ে কথা বলেছি — একটি মেথড কত কাছাকাছি প্রকৃত সমাধানে পৌঁছায়। কিন্তু কিছু সমীকরণে একটি সম্পূর্ণ ভিন্ন সমস্যা দেখা দেয়: স্টিফনেস (stiffness)। একটি সমীকরণকে "স্টিফ" বলা হয় যখন সমাধানের কোনো অংশ (বা এর একটি উপাদান) অত্যন্ত দ্রুত ক্ষয়প্রাপ্ত হয় (বড় নেতিবাচক গুণাংক), যদিও সামগ্রিক সমাধান নিজে সহজ ও মসৃণ হতে পারে। এই দ্রুত-পরিবর্তনশীল উপাদানটি ট্র্যাক করার জন্য এক্সপ্লিসিট মেথডকে অত্যন্ত ছোট step size ব্যবহার করতে বাধ্য হতে হয় — শুধু নির্ভুলতা বজায় রাখতে নয়, বরং মেথডটিকে একেবারেই ভেঙে না পড়ার জন্য।

একটি ক্লাসিক টেস্ট সমীকরণ দিয়ে এটি বিশ্লেষণ করা হয়:

এর প্রকৃত সমাধান y(x) = e^{\lambda x} — λ যত বেশি ঋণাত্মক (উদাহরণ: λ = -50), সমাধান তত দ্রুত শূন্যের দিকে যায়। এই পাঠে আমরা λ = -50 ব্যবহার করব টেস্ট সমীকরণ হিসেবে: dy/dx = -50y, y(0) = 1, y(x) = e⁻⁵⁰ˣ।

২ · অয়লার মেথডের স্ট্যাবিলিটি শর্ত

অয়লার মেথড এই টেস্ট সমীকরণে প্রয়োগ করলে:

অর্থাৎ প্রতিটি ধাপে y-কে একটি নির্দিষ্ট অ্যাম্প্লিফিকেশন ফ্যাক্টর r = 1 + hλ দিয়ে গুণ করা হয়। যদি |r| < 1, প্রতিটি ধাপে মান ছোট হতে থাকবে (স্থিতিশীল, প্রকৃত ক্ষয়ের সাথে গুণগতভাবে মিলে যায়)। কিন্তু যদি |r| > 1, প্রতিটি ধাপে মান বাড়তে থাকবে — সম্পূর্ণ অস্থির, প্রকৃত সমাধানের বিপরীত আচরণ। λ = -50 -এর জন্য স্ট্যাবিলিটি শর্ত:

অর্থাৎ h ≥ 0.04 হলে অয়লার মেথড এই সমীকরণে অস্থির হয়ে যাবে — নির্ভুলতা নয়, একেবারে স্ট্যাবিলিটিই ভেঙে পড়বে।

৩ · সত্যিকারের ডেমো — অয়লার মেথড সত্যিই বিস্ফোরিত হচ্ছে

নিচের কোড সেলে dy/dx = -50y, y(0) = 1 সমাধান করা হয়েছে h = 0.1 দিয়ে (যা স্থিতিশীলতার সীমা 0.04-এর চেয়ে অনেক বড়)। যেহেতু ইচ্ছাকৃতভাবে একটি অস্থির ডেমো দেখানো হচ্ছে, লুপে একটি ম্যাগনিটিউড গার্ড রাখা হয়েছে — মান নির্দিষ্ট সীমা পার হলে লুপ নিরাপদে থেমে যাবে, যাতে সেল কখনো hang বা crash না করে।

LAM = -50.0

def f_stiff(x, y):
    return LAM * y

def euler_stiff(h, x0=0.0, y0=1.0, xend=1.0, cap=1e8, max_steps=200):
    xs, ys = [x0], [y0]
    x, y = x0, y0
    n = min(round((xend - x0) / h), max_steps)
    blew_up = False
    for i in range(n):
        y = y + h * f_stiff(x, y)
        x = x + h
        xs.append(x); ys.append(y)
        if abs(y) > cap:      # নিরাপত্তা গার্ড -- মান অসীমের দিকে গেলে থামাও
            blew_up = True
            break
    return xs, ys, blew_up

print(f"স্ট্যাবিলিটি সীমা: h < 2/|লামডা| = 2/50 = 0.04")
print()
print("h = 0.1 (সীমার চেয়ে বড় -- অস্থির হওয়ার কথা):")
print(f"{'x':>6} | {'y':>16}")
xs, ys, blew = euler_stiff(0.1, xend=1.0)
for x, y in zip(xs, ys):
    print(f"{x:>6.2f} | {y:>16.6e}")
print("লুপ গার্ড দিয়ে থেমেছে (blew_up):", blew)

৪ · সীমানার কাছাকাছি — দোলন কিন্তু বিস্ফোরণ নয়

ঠিক সীমানায় (h = 0.04, যেখানে r = 1 + 0.04×(-50) = -1) কী হয় তা দেখা শিক্ষণীয় — মান আর বাড়েও না, কমেও না, শুধু সাইন পাল্টাতে থাকে।

print("h = 0.04 (ঠিক স্ট্যাবিলিটি সীমানায়, r = -1):")
xs, ys, blew = euler_stiff(0.04, xend=0.4)
for x, y in zip(xs, ys):
    print(f"x = {x:.3f}   y = {y:.6f}")

print()
print("h = 0.02 (সীমার ভেতরে, r = 1 + 0.02*(-50) = 0):")
xs, ys, blew = euler_stiff(0.02, xend=0.1)
for x, y in zip(xs, ys):
    print(f"x = {x:.3f}   y = {y:.6f}")

৫ · একটি সত্যিকারের স্থিতিশীল সমাধান

এবার স্থিতিশীলতা সীমার ভেতরে একটি বাস্তবসম্মত h ব্যবহার করে দেখা যাক মেথডটি সত্যিই কেমন কাজ করে — প্রকৃত ক্লোজড-ফর্ম সমাধান e⁻⁵⁰ˣ-এর বিপরীতে।

import math

def exact_stiff(x):
    return math.exp(LAM * x)

print("h = 0.01 (স্থিতিশীল সীমার অনেক ভেতরে):")
print(f"{'x':>6} | {'y_euler':>14} | {'y_exact':>14}")
xs, ys, blew = euler_stiff(0.01, xend=0.3)
for x, y in zip(xs, ys):
    print(f"{x:>6.2f} | {y:>14.6e} | {exact_stiff(x):>14.6e}")

final_err = abs(ys[-1] - exact_stiff(xs[-1]))
print()
print(f"চূড়ান্ত বিন্দুতে প্রকৃত এরর: {final_err:.3e}")

৬ · RK4-ও অস্থির হতে পারে, তবে বড় সীমায়

এটি একটি সাধারণ ভুল ধারণা যে "উচ্চতর-ক্রম মেথড মানেই স্ট্যাবিলিটির সমস্যা নেই" — RK4-এরও একটি সসীম স্ট্যাবিলিটি অঞ্চল আছে, শুধু তা অয়লারের চেয়ে বড়। নিচের কোড সেলে RK4 একই স্টিফ সমীকরণে h = 0.1 -এ চালিয়ে দেখানো হয়েছে এটিও অস্থির হয়ে যায়, এবং তারপর কোন h-এ RK4 স্থিতিশীল থাকে তা পরীক্ষা করা হয়েছে।

def rk4_stiff(h, x0=0.0, y0=1.0, xend=1.0, cap=1e6, max_steps=60):
    x, y = x0, y0
    n = min(round((xend - x0) / h), max_steps)
    blew_up = False
    for i in range(n):
        k1 = f_stiff(x, y)
        k2 = f_stiff(x + h / 2, y + h / 2 * k1)
        k3 = f_stiff(x + h / 2, y + h / 2 * k2)
        k4 = f_stiff(x + h, y + h * k3)
        y = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
        x = x + h
        if abs(y) > cap:
            blew_up = True
            break
    return y, blew_up, i + 1

print("RK4, h = 0.1 (অয়লারের সীমা 0.04-এর চেয়ে বড়, RK4-এরও সীমা পার হতে পারে):")
yf, blew, steps = rk4_stiff(0.1, xend=1.0)
print(f"শেষ y = {yf:.4e}   blew_up = {blew}   ধাপ সম্পন্ন = {steps}")
print()

print("বিভিন্ন h-এ RK4 স্থিতিশীল থাকে কিনা:")
for h in [0.04, 0.05, 0.06, 0.07, 0.08]:
    yf, blew, steps = rk4_stiff(h, xend=0.5, cap=1e6, max_steps=50)
    print(f"h = {h:.2f}   শেষ y = {yf:>14.6e}   blew_up = {blew}")