পদার্থবিজ্ঞান ও ইঞ্জিনিয়ারিং সিমুলেশনে নিউমেরিক্যাল মেথডস

১ · সিমুলেশন মানেই আসলে ODE সলভিং

পদার্থবিজ্ঞান ও ইঞ্জিনিয়ারিং-এ যখন আমরা "সিমুলেট করা" বলি — একটি সেতুর কম্পন, একটি গাড়ির সাসপেনশন সিস্টেম, একটি স্যাটেলাইটের কক্ষপথ, একটি রোবট আর্মের গতি — প্রায় প্রতিটি ক্ষেত্রেই আসল গাণিতিক সমস্যাটি একটি ডিফারেনশিয়াল সমীকরণ যা সিস্টেমের অবস্থা (position, velocity, temperature ইত্যাদি) সময়ের সাথে কীভাবে পরিবর্তিত হয় তা বর্ণনা করে। বেশিরভাগ বাস্তব সিস্টেমের জন্য এই সমীকরণের কোনো ক্লোজড-ফর্ম বীজগাণিতিক সমাধান নেই (M6-এর ভূমিকায় যেমন বলা হয়েছিল) — তাই M6-M7-এ শেখা Euler, Heun ও RK4-এর মতো নিউমেরিক্যাল ODE সলভারই এই ক্ষেত্রের আসল কর্মী হাতিয়ার।

২ · স্প্রিং-মাস সিস্টেম — গতির সমীকরণ

একটি ভর m স্প্রিং (স্প্রিং ধ্রুবক k) ও ড্যাম্পার (ড্যাম্পিং সহগ c, যেমন ঘর্ষণ বা বাতাসের রোধ) দিয়ে সংযুক্ত থাকলে, নিউটনের দ্বিতীয় সূত্র থেকে সিস্টেমের গতির সমীকরণ:

$$ m\,x''(t) + c\,x'(t) + k\,x(t) = 0 $$

এটি একটি দ্বিতীয়-ক্রম ODE — কিন্তু M6/L27-এ শেখা কৌশল ব্যবহার করে একে দুটি প্রথম-ক্রম ODE-এর একটি সিস্টেমে রূপান্তর করা যায়, স্টেট ভেক্টর y = [position, velocity] ব্যবহার করে:

$$ x' = v, \qquad v' = \frac{-k\,x - c\,v}{m} $$

এই দুটি সমীকরণ একসাথে RK4-এর মাধ্যমে সময়ের প্রতিটি ছোট ধাপে সমাধান করা যায় — ঠিক M6/L30-এ শেখা RK4 প্যাটার্নটিই এখানে দুই-উপাদান স্টেট ভেক্টরে প্রয়োগ করা হয়েছে।

৩ · সত্যিকারের সিমুলেশন — RK4 দিয়ে ড্যাম্পড দোলন

নিচের কোড সেলে m = 1.0 kg, k = 4.0 N/m, ড্যাম্পিং সহগ c = 0.4, প্রাথমিক পজিশন x₀ = 1.0 m ও প্রাথমিক ভেলোসিটি v₀ = 0 নিয়ে, স্টেপ সাইজ h = 0.1 সেকেন্ড দিয়ে t = 0 থেকে t = 3 সেকেন্ড পর্যন্ত RK4 দিয়ে সিস্টেমটি ইন্টিগ্রেট করা হয়েছে — সম্পূর্ণ Python list-ভিত্তিক, কোনো NumPy ছাড়াই।

import math

# Spring-mass system: m*x'' + c*x' + k*x = 0
# State vector y = [position, velocity]
m = 1.0      # kg
k = 4.0      # N/m
c = 0.4      # damping coefficient (kg/s)

def derivatives(t, y):
    x, v = y
    a = (-k * x - c * v) / m
    return [v, a]

def rk4_step(t, y, h):
    k1 = derivatives(t, y)
    y2 = [y[i] + h / 2 * k1[i] for i in range(2)]
    k2 = derivatives(t + h / 2, y2)
    y3 = [y[i] + h / 2 * k2[i] for i in range(2)]
    k3 = derivatives(t + h / 2, y3)
    y4 = [y[i] + h * k3[i] for i in range(2)]
    k4 = derivatives(t + h, y4)
    return [y[i] + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) for i in range(2)]

x0, v0 = 1.0, 0.0
h = 0.1
T = 3.0
n_steps = int(T / h)

t = 0.0
y = [x0, v0]

print(f"{'t':>5} | {'x (position)':>14} | {'v (velocity)':>14}")
for i in range(n_steps + 1):
    if i % 5 == 0:
        print(f"{t:5.2f} | {y[0]:14.6f} | {y[1]:14.6f}")
    if i < n_steps:
        y = rk4_step(t, y, h)
        t += h

print()
print(f"চূড়ান্ত অবস্থা t={t:.2f}: x={y[0]:.6f}  v={y[1]:.6f}")

৪ · সত্যতা যাচাই — আনড্যাম্পড কেসে ক্লোজড-ফর্ম সমাধানের সাথে তুলনা

একটি সিমুলেশন "কাজ করছে বলে মনে হওয়া" আর "সত্যিই সঠিক হওয়া" — দুটো ভিন্ন জিনিস। যেকোনো নতুন সিমুলেশন কোডকে বিশ্বাস করার আগে অন্তত একটি এমন বিশেষ ক্ষেত্রে যাচাই করা উচিত যেখানে সঠিক উত্তর আগে থেকেই জানা আছে — এটাই M11-এর "যাচাইযোগ্যতা" নীতির সরাসরি প্রয়োগ। যখন ড্যাম্পিং c = 0, স্প্রিং-মাস সিস্টেমের একটি পরিচিত ক্লোজড-ফর্ম সমাধান আছে:

$$ x(t) = x_0 \cos(\omega t), \qquad \omega = \sqrt{k/m} $$

নিচের কোড সেলে একই RK4 সিমুলেশন c = 0 দিয়ে চালিয়ে, প্রতিটি টাইম-স্টেপে RK4-এর ফলাফল এই অ্যানালিটিক সূত্রের সাথে তুলনা করা হয়েছে — সত্যিকারের কম্পিউটেড এরর দেখতে।

import math

m, k = 1.0, 4.0
c2 = 0.0   # undamped case

def derivatives_undamped(t, y):
    x, v = y
    a = (-k * x - c2 * v) / m
    return [v, a]

def rk4_step_undamped(t, y, h):
    k1 = derivatives_undamped(t, y)
    y2 = [y[i] + h / 2 * k1[i] for i in range(2)]
    k2 = derivatives_undamped(t + h / 2, y2)
    y3 = [y[i] + h / 2 * k2[i] for i in range(2)]
    k3 = derivatives_undamped(t + h / 2, y3)
    y4 = [y[i] + h * k3[i] for i in range(2)]
    k4 = derivatives_undamped(t + h, y4)
    return [y[i] + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) for i in range(2)]

x0, v0 = 1.0, 0.0
h = 0.1
n_steps = 30
omega = math.sqrt(k / m)
print(f"omega = sqrt(k/m) = {omega:.6f}")
print()

t = 0.0
y = [x0, v0]
print(f"{'t':>5} | {'RK4 x':>12} | {'analytic x':>12} | {'abs error':>12}")
for i in range(n_steps + 1):
    analytic = x0 * math.cos(omega * t)
    err = abs(y[0] - analytic)
    if i % 5 == 0:
        print(f"{t:5.2f} | {y[0]:12.8f} | {analytic:12.8f} | {err:12.2e}")
    if i < n_steps:
        y = rk4_step_undamped(t, y, h)
        t += h

final_analytic = x0 * math.cos(omega * t)
final_err = abs(y[0] - final_analytic)
print()
print(f"চূড়ান্ত RK4 x = {y[0]:.8f}, অ্যানালিটিক x = {final_analytic:.8f}, প্রকৃত এরর = {final_err:.2e}")