একটি নিউমেরিক্যাল সলভার টুলকিট তৈরি — ডিজাইন
এই পাঠে যা শিখবেন
- কীভাবে একটি নির্দিষ্ট সমীকরণে বাঁধা কোড থেকে একটি জেনেরিক, পুনর্ব্যবহারযোগ্য ফাংশন ইন্টারফেসে যাওয়া
যায় (যেমন
rk4_step(f, t, y, h)— যেকোনোf-এর জন্য কাজ করে) - M11/L52-এর সিদ্ধান্ত-কাঠামো বাস্তবে প্রয়োগ করে একটি টুলকিটের প্রতিটি ফাংশনের জন্য মেথড বাছাই করা
- একটি ODE সলভার ও একটি রুট-ফাইন্ডারকে একসাথে চেইন করে একটি জটিলতর প্রশ্নের উত্তর বের করা
- L57-এর ক্যাপস্টোনের জন্য সমস্যা ও টুলকিট প্রস্তুত করা
১ · ক্যাপস্টোন সমস্যা বাছাই — নিউটনের শীতলীকরণ সূত্র
একটি টুলকিট ডিজাইন করার আগে জানা দরকার এটি ঠিক কী সমস্যা সমাধান করবে। এই মডিউলের জন্য আমরা একটি
সম্পূর্ণ হাইপোথেটিক্যাল ও ইলাস্ট্রেটিভ ইঞ্জিনিয়ারিং সিনারিও বেছে নিচ্ছি — বাস্তব পরিমাপ করা
ডেটা নয়, একটি teaching example: একটি ধাতব যন্ত্রাংশ একটি ফার্নেস থেকে T₀ = 180°C তাপমাত্রায়
বের করা হয়েছে এবং একটি ঘরে রাখা হয়েছে যার তাপমাত্রা T_env = 22°C। প্রকৌশলীরা জানতে চান পরবর্তী
অপারেশন (যেমন কোটিং) শিডিউল করার আগে যন্ত্রাংশটি নিরাপদে হাত দেওয়ার তাপমাত্রা T_safe = 45°C-এ
পৌঁছাতে কত সময় লাগবে।
নিউটনের শীতলীকরণ সূত্র অনুযায়ী, তাপমাত্রা পরিবর্তনের হার তাপমাত্রা ও পরিবেশের তাপমাত্রার পার্থক্যের সমানুপাতিক:
$$\frac{dT}{dt} = -k(T - T_{env})$$
এখানে k একটি শীতলীকরণ ধ্রুবক (উপাদান, পৃষ্ঠতল ও বাতাস চলাচলের উপর নির্ভর করে) — আমরা এই মডিউলে
একটি ইলাস্ট্রেটিভ মান k = 0.045 প্রতি মিনিট ব্যবহার করব। লক্ষ্য করুন এটি M6/M7-এর একটি সাধারণ
প্রাথমিক-মূল্য সমস্যা (IVP) — এবং এর একটি ক্লোজড-ফর্ম সমাধানও আছে (T(t) = T_env + (T₀-T_env)e^{-kt}),
যা L57-এ ভ্যালিডেশনের জন্য কাজে আসবে। কিন্তু এখানে আমরা ইচ্ছাকৃতভাবে একে নিউমেরিক্যালি সমাধান
করব, কারণ একটি টুলকিট এমনভাবে ডিজাইন করা উচিত যা এমন সমীকরণেও কাজ করবে যাদের কোনো ক্লোজড-ফর্ম সমাধান নেই।
২ · টুলকিট ডিজাইন — কোন ফাংশন, কেন
একটি ভালো টুলকিট ডিজাইনের প্রথম প্রশ্ন: কোন কোন বিল্ডিং-ব্লক লাগবে? এই সমস্যার জন্য আমরা তিনটি ফাংশন-মডিউল চিহ্নিত করেছি:
rk4_step(f, t, y, h) ও solve_ode(f, t0, y0, h, n) — যেকোনো dy/dt = f(t,y) আকারের সমীকরণ ধাপে ধাপে সমাধান করে, শুধু শীতলীকরণ সমীকরণে বাঁধা নয়।bisection_root(g, a, b, tol) — যেকোনো g(x)=0 সমাধান করে; আমরা এটি g(t) = ODE_solution(t) - T_safe-এ প্রয়োগ করব, অর্থাৎ ODE সলভারের আউটপুটের উপর রুট-ফাইন্ডিং চালাব।একটি ছোট হেল্পার যা একই সমস্যা দুটি ভিন্ন স্টেপ-সাইজে সমাধান করে এবং এরর-রেশিও গণনা করে — M5/L22-এর রিচার্ডসন-ধাঁচের ধারণা, L57-এ কনভারজেন্স-অর্ডার যাচাই করতে ব্যবহৃত হবে।
M11/L52-এর সিদ্ধান্ত-কাঠামো অনুযায়ী, ODE সলভার বাছাইয়ের মূল প্রশ্ন হলো অ্যাকুরেসি বনাম কম্পিউটেশনাল খরচের
ট্রেড-অফ। M6/L28-এ দেখা গিয়েছিল Euler মেথড প্রথম-অর্ডার (এরর O(h)) এবং M6/L30-এ RK4
চতুর্থ-অর্ডার (এরর O(h⁴))। যেহেতু আমাদের শীতলীকরণ সমীকরণটি "smooth" (কোনো হঠাৎ পরিবর্তন বা
stiffness নেই — M6/L32-এর অর্থে), RK4 প্রতি ধাপে সামান্য বেশি গণনা করেও অনেক বড় স্টেপ-সাইজে একই অ্যাকুরেসি
দেয় — তাই টুলকিটের ডিফল্ট ODE সলভার হিসেবে RK4 বাছাই করা হয়েছে।
M2/L09-এর তুলনা-টেবিল অনুযায়ী নিউটন-রাফসন (M2/L07) সাধারণত বাইসেকশনের চেয়ে দ্রুত কনভার্জ করে (কোয়াড্রাটিক
বনাম লিনিয়ার) — কিন্তু নিউটন-রাফসনের জন্য g'(x) (ডেরিভেটিভ) দরকার। আমাদের
g(t) = solve_to_time(t) - T_safe ফাংশনটি একটি সম্পূর্ণ RK4 সিমুলেশনের আউটপুট — এর কোনো
সহজ, বিশ্লেষণাত্মক ডেরিভেটিভ নেই (এটিকে সংখ্যাগতভাবে অনুমান করা যেত, কিন্তু তাতে জটিলতা ও অতিরিক্ত এরর যোগ
হতো)। M11/L52-এর কাঠামো অনুযায়ী: যখন ডেরিভেটিভ সহজলভ্য নয় কিন্তু ফাংশনটি নির্ভরযোগ্যভাবে সাইন-পরিবর্তন করে,
বাইসেকশনের নিশ্চয়তা ও দৃঢ়তা তার ধীরগতির চেয়ে বেশি গুরুত্বপূর্ণ — তাই টুলকিটের ডিফল্ট
রুট-ফাইন্ডার হিসেবে বাইসেকশন বাছাই করা হয়েছে।
৩ · কোড সেল — টুলকিট ফাংশন লেখা ও টেস্ট করা
নিচের কোড সেলে দুটি জেনেরিক ফাংশন লেখা ও টেস্ট করা হয়েছে: প্রথমে rk4_step/solve_ode
শীতলীকরণ সমীকরণে প্রয়োগ করে একটি তাপমাত্রা-ট্রেস দেখানো হয়েছে; তারপর bisection_root একটি সাধারণ
টেস্ট-ফাংশনে যাচাই করা হয়েছে; এবং শেষে দুটো ফাংশন একসাথে চেইন করে t_safe (৪৫°C-এ পৌঁছাতে সময়)
বের করা হয়েছে — এই চেইনিং প্যাটার্নটিই L57-এর ক্যাপস্টোনের মূল কৌশল।
import math
# ---------- জেনেরিক RK4 স্টেপার (M6/L30-এর ধারণা, জেনেরিক ইন্টারফেসে) ----------
def rk4_step(f, t, y, h):
k1 = f(t, y)
k2 = f(t + h/2, y + h/2 * k1)
k3 = f(t + h/2, y + h/2 * k2)
k4 = f(t + h, y + h * k3)
return y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
def solve_ode(f, t0, y0, h, n_steps):
t, y = t0, y0
trace = [(t, y)]
for _ in range(n_steps):
y = rk4_step(f, t, y, h)
t = t + h
trace.append((t, y))
return trace
# টেস্ট: নিউটনের শীতলীকরণ সূত্র dT/dt = -k(T - T_env)
T_ENV = 22.0
K = 0.045
def cooling_ode(t, T):
return -K * (T - T_ENV)
trace = solve_ode(cooling_ode, 0.0, 180.0, 5.0, 12)
for t, T in trace:
print(f"t={t:5.1f} মিনিট T={T:7.3f}°C")
# ---------- জেনেরিক বাইসেকশন রুট-ফাইন্ডার (M2/L05-এর ধারণা, জেনেরিক ইন্টারফেসে) ----------
def bisection_root(g, a, b, tol=1e-6, max_iter=100):
ga, gb = g(a), g(b)
if ga * gb > 0:
raise ValueError("no sign change in [a,b]")
for i in range(max_iter):
mid = (a + b) / 2
gm = g(mid)
if abs(gm) < tol or (b - a) / 2 < tol:
return mid, i + 1
if ga * gm < 0:
b = mid
else:
a, ga = mid, gm
return (a + b) / 2, max_iter
# স্যানিটি টেস্ট -- একটি সাধারণ ফাংশনে
def test_f(x):
return x**3 - x - 2
root, iters = bisection_root(test_f, 1.0, 2.0, tol=1e-8)
print(f"\nটেস্ট রুট (x^3-x-2): {root:.8f}, {iters} ইটারেশনে, f(root)={test_f(root):.2e}")
# টুলকিট চেইনিং: বাইসেকশন + RK4 -- ঠান্ডা হয়ে ৪৫°C-এ পৌঁছাতে কত সময় লাগে?
def solve_to_time(t_target, h):
n_steps = max(1, round(t_target / h))
h_actual = t_target / n_steps
t, y = 0.0, 180.0
for _ in range(n_steps):
y = rk4_step(cooling_ode, t, y, h_actual)
t += h_actual
return y
def g_target(t):
return solve_to_time(t, 1.0) - 45.0
t_root, iters2 = bisection_root(g_target, 0.0, 120.0, tol=1e-6)
print(f"\n৪৫°C-এ পৌঁছাতে সময় (বাইসেকশন+RK4): t={t_root:.6f} মিনিট, {iters2} ইটারেশনে")
print(f"যাচাই: t={t_root:.6f}-এ T = {solve_to_time(t_root,1.0):.6f}°C")
t=0-এ ১৮০°C থেকে দ্রুত কমে
t=60 মিনিটে ৩২.৬১৯°C-এ পৌঁছেছে। টেস্ট-ফাংশনের রুট ১.৫২১৩৭৯৭০
মাত্র ২৭ ইটারেশনে পাওয়া গেছে (f(root) ≈ -2.98×10⁻⁸, কার্যত শূন্য)। সবচেয়ে
গুরুত্বপূর্ণ: দুটো ফাংশন চেইন করে বের করা t_safe = 42.824464 মিনিট — অর্থাৎ যন্ত্রাংশটি
নিরাপদে হাত দেওয়ার তাপমাত্রায় পৌঁছাতে প্রায় ৪২.৮ মিনিট লাগবে, এবং যাচাই-প্রিন্ট নিশ্চিত করছে
সেই সময়ে প্রকৃতপক্ষে তাপমাত্রা ঠিক ৪৫.০০০০০০°C।
৪ · ডিজাইন সিদ্ধান্ত ও ট্রেড-অফের সারসংক্ষেপ
লক্ষ্য করুন আমরা solve_to_time-কে ইচ্ছাকৃতভাবে "সম্পূর্ণ ODE সলভ করে একটি নির্দিষ্ট সময়ের
তাপমাত্রা ফেরত দেয়" — এই আকারে ডিজাইন করেছি, যাতে এটি সরাসরি bisection_root-এর
g(x) আর্গুমেন্ট হিসেবে বসানো যায়। এটাই একটি ভালো টুলকিট ডিজাইনের মূল নীতি: প্রতিটি ফাংশনের
ইনপুট/আউটপুট এমনভাবে ডিজাইন করা যাতে সেগুলো একে অপরের সাথে সহজে চেইন করা যায়, প্রতিবার নতুন
করে "গ্লু কোড" না লিখে। L57-এর ক্যাপস্টোনে আমরা ঠিক এই তিনটি বিল্ডিং-ব্লক (ODE সলভার, রুট-ফাইন্ডার, ও
এরর-চেকার) ব্যবহার করব, তবে আরও কঠোরভাবে — একটি ক্লোজড-ফর্ম রেফারেন্সের বিপরীতে ভ্যালিডেট করে, স্টেপ-সাইজ
অর্ধেক করে কনভারজেন্স-অর্ডার সত্যিই যাচাই করে, এবং k-এর অনিশ্চয়তার প্রভাব একটি মন্টে কার্লো
সেনসিটিভিটি অ্যানালাইসিস (M11/L51-এর প্যাটার্নে) দিয়ে পরিমাপ করে।
একটি নিউমেরিক্যাল টুলকিট ডিজাইন করা মানে শুধু "কোন মেথড কাজ করে" তা জানা নয় — এটি এমন ফাংশন-ইন্টারফেস ডিজাইন করা যা জেনেরিক (একটি নির্দিষ্ট সমস্যায় বাঁধা নয়) এবং চেইনযোগ্য (একটির আউটপুট আরেকটির ইনপুট হতে পারে)। প্রতিটি মেথড বাছাই M11/L52-এর সিদ্ধান্ত-কাঠামো অনুযায়ী যুক্তিসহ নেওয়া উচিত — অ্যাকুরেসি, দৃঢ়তা, ও কম্পিউটেশনাল খরচের ট্রেড-অফ বিবেচনা করে, শুধু "যা পরিচিত তাই ব্যবহার করা" নয়।
ভাবনার প্রশ্ন
প্রতিটি প্রশ্ন নিজে কিছুক্ষণ ভাবুন — তারপর "→ উত্তর" চাপুন।
প্র ০১
solve_to_time ফাংশনটি প্রতিবার t=0 থেকে সম্পূর্ণ ODE পুনরায় সমাধান করে,
আগের ফলাফল পুনর্ব্যবহার করে না। বাইসেকশনের ২৩টি ইটারেশনের প্রতিটিতে এটি ঘটে। এটি কি একটি ডিজাইন সমস্যা?
হ্যাঁ, কম্পিউটেশনাল দৃষ্টিকোণ থেকে এটি অদক্ষ — একই সমীকরণ ২৩ বার সমাধান করা হচ্ছে যখন একটি স্মার্ট ডিজাইন একবার সমাধান করে ফলাফল ক্যাশ করে বাইসেকশনের প্রতিটি ধাপে ইন্টারপোলেট করতে পারত। কিন্তু এই ছোট সমস্যায় (প্রতিটি সলভ মাত্র কয়েক মিলিসেকেন্ড) সরলতা অগ্রাধিকার পেয়েছে — এটি একটি বাস্তব ট্রেড-অফ: বড় বা ব্যয়বহুল সিমুলেশনে এই ক্যাশিং অপ্টিমাইজেশন অপরিহার্য হয়ে উঠত (M10/L49-এর কম্পিউটেশনাল-খরচ আলোচনার সাথে সংযুক্ত)।
প্র ০২ যদি ভবিষ্যতে টুলকিটে একটি "স্টিফ" ODE (M6/L32) যোগ করতে হয়, RK4 কি এখনও ডিফল্ট বাছাই থাকবে?
না, সম্ভবত না। M6/L32-এ দেখা গিয়েছিল স্টিফ সমীকরণে এক্সপ্লিসিট মেথড (Euler, RK4 উভয়ই এক্সপ্লিসিট)
স্টেবিলিটি বজায় রাখতে অত্যন্ত ছোট স্টেপ-সাইজ দাবি করে, যা কম্পিউটেশনালি ব্যয়বহুল হয়ে যায়। একটি
দৃঢ় টুলকিট ডিজাইনে solve_ode-এর পাশে একটি ইমপ্লিসিট মেথডও (M6/L32-এ উল্লেখিত ধারণা)
বিকল্প হিসেবে রাখা উচিত, এবং সমস্যার প্রকৃতি অনুযায়ী কোনটি ব্যবহার করতে হবে তা M11/L52-এর কাঠামো দিয়ে
সিদ্ধান্ত নেওয়া উচিত।
প্র ০৩
bisection_root ফাংশনটি ga * gb > 0 হলে একটি এরর তোলে (কোনো সাইন-চেঞ্জ
নেই)। কেন এই চেক টুলকিট ডিজাইনের একটি গুরুত্বপূর্ণ অংশ?
কারণ এটি "সাইলেন্ট ফেইলিওর" প্রতিরোধ করে — L01-এর ভাবনার প্রশ্নেও এই একই নীতি এসেছিল। একটি পুনর্ব্যবহারযোগ্য টুলকিট ফাংশন যদি ভুল ইনপুটে চুপচাপ ভুল উত্তর দেয়, তাহলে পরবর্তী কোনো লেসনে (যেমন L57-এ) সেই ফাংশনটি ভুলভাবে ব্যবহার করলে ভুল ধরা কঠিন হয়ে যায়। একটি ভালো টুলকিট ডিজাইনের নীতি হলো preconditions স্পষ্টভাবে চেক করে ব্যর্থ হওয়া (fail loudly), অনুমান করে চুপচাপ এগিয়ে যাওয়া নয়।
অনুশীলন
-
চিন্তা করুন: কোড সেলে
bisection_root(g_target, 0.0, 120.0, tol=1e-6)-এরtolমান1e-6-এর বদলে1e-10করলে ইটারেশন সংখ্যা মোটামুটি কতটা বাড়বে বলে আপনার ধারণা?বাইসেকশন প্রতি ইটারেশনে ইন্টারভাল অর্ধেক করে (L01/M2-এর লিনিয়ার কনভারজেন্স), তাই টলারেন্স
10⁴গুণ কমাতে (1e-6 → 1e-10) আনুমানিকlog₂(10⁴) ≈ 13.3টি অতিরিক্ত ইটারেশন লাগবে — অর্থাৎ মোট ইটারেশন প্রায়23 + 13 ≈ 36-এর কাছাকাছি হওয়া উচিত। -
পরীক্ষা করুন: কোড সেলে
tol=1e-6-কেtol=1e-10-এ পরিবর্তন করে Run চেপে আপনার অনুমান যাচাই করুন।ফলাফল নিশ্চিত করবে ইটারেশন সংখ্যা বেড়ে প্রায় ৩৬-৩৭-এ পৌঁছায় (ঠিক সংখ্যাটি ইন্টারভাল প্রস্থের উপর নির্ভর করে সামান্য ভিন্ন হতে পারে), এবং
t_root-এর মান আরও কয়েক দশমিক স্থান পর্যন্ত সঠিক হয় — কিন্তু ৪২.৮২ মিনিটের কাছাকাছি মূল উত্তর কার্যত অপরিবর্তিত থাকে, কারণtol=1e-6নিজেই ইতিমধ্যে অনেক বেশি নির্ভুল ছিল এই বাস্তব সমস্যার জন্য।
আরও পড়ুন · ABCL TECH-এ আপনার পরবর্তী পদক্ষেপ
- ক্যাপস্টোন — এই টুলকিট দিয়ে সমস্যাটি এন্ড-টু-এন্ড সমাধান করুন পরবর্তী পাঠ এই পাঠে ডিজাইন করা টুলকিট ব্যবহার করে ভ্যালিডেশন, কনভারজেন্স-অর্ডার চেক, ও মন্টে কার্লো সেনসিটিভিটি অ্যানালাইসিস সহ সম্পূর্ণ সমস্যাটি সমাধান করা হবে।
- M11/L52 — সঠিক মেথড বাছাই: সিদ্ধান্ত-কাঠামো রিভিশন এই পাঠের দুটি ডিজাইন-সিদ্ধান্তই যে কাঠামো থেকে এসেছে, তা আরও বিস্তারিতভাবে দেখুন।
- কোর্সের সম্পূর্ণ সিলেবাস দেখুন ৫৭টি পাঠ এরর অ্যানালাইসিস থেকে ক্যাপস্টোন পর্যন্ত — সম্পূর্ণ কোর্স ম্যাপ।