মাল্টিস্টেপ মেথড — অ্যাডামস-ব্যাশফোর্থ/মৌলটন

১ · সিঙ্গল-স্টেপ বনাম মাল্টিস্টেপ

L27–L30-এর প্রতিটি মেথড — অয়লার, হয়েন, RK4 — সিঙ্গল-স্টেপ: y_{n+1} বের করতে শুধু (x_n, y_n)-এর তথ্য লাগে, তার আগের কোনো বিন্দুর তথ্য নয়। RK4 প্রতি ধাপে চারবার f মূল্যায়ন করে নির্ভুলতা কেনে — কিন্তু প্রতিটি মূল্যায়নই "একবার ব্যবহার করে ফেলে দেওয়া" হয়।

মাল্টিস্টেপ মেথড একটি ভিন্ন কৌশল নেয়: যেহেতু আমরা ইতিমধ্যে আগের কয়েকটি বিন্দুতে f-এর মান গণনা করে ফেলেছি (আগের ধাপগুলো থেকে), কেন সেই তথ্য ফেলে দেওয়া হবে? এর বদলে, আগের কয়েকটি বিন্দুর f-মান দিয়ে একটি পলিনোমিয়াল ফিট করে (M4-এর ইন্টারপোলেশনের সাথে ধারণাগতভাবে ঘনিষ্ঠ) সেই পলিনোমিয়াল ব্যবহার করে পরবর্তী বিন্দু আনুমানিক করা হয় — নতুন করে অতিরিক্ত মূল্যায়ন ছাড়াই (শুধু বর্তমান বিন্দুর একটি নতুন মূল্যায়ন যোগ হয়)।

২ · ২-স্টেপ অ্যাডামস-ব্যাশফোর্থ (AB2)

সবচেয়ে সরল এক্সপ্লিসিট মাল্টিস্টেপ মেথডগুলোর একটি হলো ২-স্টেপ অ্যাডামস-ব্যাশফোর্থ, যা শেষ দুই বিন্দুর f-মান দিয়ে একটি রৈখিক এক্সট্রাপোলেশন ব্যবহার করে:

লক্ষ্য করুন — এই সূত্রে মাত্র একটি নতুন f-মূল্যায়ন লাগে (f(x_n, y_n)); f(x_{n-1}, y_{n-1}) আগের ধাপ থেকে ইতিমধ্যে গণনা করা এবং সংরক্ষিত আছে। ওজন 3/2 ও -1/2 এমনভাবে নির্বাচিত (রৈখিক এক্সট্রাপোলেশনের গাণিতিক ফলাফল) যাতে সামগ্রিক মেথড দ্বিতীয়-ক্রম নির্ভুল (O(h²)) হয় — অর্থাৎ প্রতি ধাপে RK4-এর চতুর্থাংশ হিসাব-খরচে (১টি বনাম ৪টি f-মূল্যায়ন) হয়েনের সমান কনভারজেন্স-অর্ডার পাওয়া যায়।

সমস্যা: y_{n+1} বের করতে y_n ও y_{n-1} দুটোই লাগে — কিন্তু শুরুতে আমাদের কাছে শুধু y_0 (শুরুর শর্ত) আছে, y_1 নেই। তাই প্রথম ধাপটি এই মেথড নিজে চালাতে পারে না — একে বলে বুটস্ট্র্যাপিং সমস্যা, এবং এর সমাধান হলো প্রথম ধাপের জন্য একটি নির্ভরযোগ্য সিঙ্গল-স্টেপ মেথড (RK4) ব্যবহার করা, তারপর থেকে AB2 নিজেই চালিয়ে যাওয়া।

৩ · সত্যিকারের বাস্তবায়ন — RK4 দিয়ে বুটস্ট্র্যাপ করে AB2

L27–L30-এর একই টেস্ট IVP (dy/dx = -y, y(0) = 1, y(x) = e⁻ˣ) ব্যবহার করে, নিচের কোড সেলে AB2 হাতে-লেখা হয়েছে — প্রথম ধাপ RK4 দিয়ে, তারপরের সব ধাপ AB2-এর নিজস্ব সূত্র দিয়ে।

import math

def f(x, y):
    return -y

def exact(x):
    return math.exp(-x)

def rk4_step(x, y, h):
    k1 = f(x, y)
    k2 = f(x + h / 2, y + h / 2 * k1)
    k3 = f(x + h / 2, y + h / 2 * k2)
    k4 = f(x + h, y + h * k3)
    return y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)

def ab2(h, x0=0.0, y0=1.0, xend=2.0):
    xs, ys = [x0], [y0]
    n = round((xend - x0) / h)

    # বুটস্ট্র্যাপ: প্রথম ধাপ RK4 দিয়ে (AB2 নিজে এটা করতে পারে না)
    x1 = x0 + h
    y1 = rk4_step(x0, y0, h)
    xs.append(x1); ys.append(y1)

    f_prev = f(x0, y0)   # f_{n-1}
    f_curr = f(x1, y1)   # f_n
    x, y = x1, y1
    for i in range(1, n):
        y_next = y + h * (1.5 * f_curr - 0.5 * f_prev)   # AB2 সূত্র
        x_next = x + h
        xs.append(x_next); ys.append(y_next)
        f_prev = f_curr
        f_curr = f(x_next, y_next)
        x, y = x_next, y_next
    return xs, ys

print("AB2 দিয়ে dy/dx = -y, y(0) = 1 সমাধান, h = 0.2:")
print(f"{'x':>6} | {'y_ab2':>10} | {'y_exact':>10} | {'|error|':>10}")
xs, ys = ab2(0.2)
for x, y in zip(xs, ys):
    err = abs(y - exact(x))
    print(f"{x:>6.2f} | {y:>10.6f} | {exact(x):>10.6f} | {err:>10.6f}")

৪ · সত্যিকারের তুলনা — AB2 বনাম RK4

নিচের কোড সেলে একই h-মানে AB2 ও RK4-এর চূড়ান্ত এরর সরাসরি পাশাপাশি তুলনা করা হয়েছে।

def rk4(h, x0=0.0, y0=1.0, xend=2.0):
    x, y = x0, y0
    n = round((xend - x0) / h)
    for i in range(n):
        y = rk4_step(x, y, h)
        x = x + h
    return y

ex = exact(2.0)
hs = [0.4, 0.2, 0.1, 0.05, 0.025]
print(f"{'h':>8} | {'err_ab2':>12} | {'err_rk4':>16} | {'ab2/rk4':>10}")
prev_ab2 = None
for h in hs:
    ab2_final = ab2(h)[1][-1]
    rk4_final = rk4(h)
    err_ab2 = abs(ab2_final - ex)
    err_rk4 = abs(rk4_final - ex)
    print(f"{h:>8} | {err_ab2:>12.8f} | {err_rk4:>16.10f} | {err_ab2 / err_rk4:>10.1f}")

print()
print("AB2-এর নিজস্ব কনভারজেন্স-অর্ডার যাচাই (h অর্ধেক -> এরর প্রায় কত ভাগ?):")
prev_err = None
for h in hs:
    err = abs(ab2(h)[1][-1] - ex)
    ratio = (prev_err / err) if prev_err is not None else None
    ratio_str = f"{ratio:.3f}" if ratio is not None else "  --  "
    print(f"h={h:<8} err={err:.8f}  ratio_prev/curr={ratio_str}")
    prev_err = err