ক্যাপস্টোন — একটি বাস্তব-বিশ্বের ইঞ্জিনিয়ারিং সমস্যা নিউমেরিক্যালি এন্ড-টু-এন্ড সমাধান

# ========== ধাপ ১: গভর্নিং ইকুয়েশন সেটআপ (হাইপোথেটিক্যাল সিনারিও) ==========
T_ENV = 22.0     # ঘরের তাপমাত্রা, °C
T0 = 180.0       # t=0-এ যন্ত্রাংশের তাপমাত্রা, °C
K = 0.045        # শীতলীকরণ ধ্রুবক, প্রতি মিনিট (ইলাস্ট্রেটিভ মান)
T_SAFE = 45.0    # নিরাপদে হাত দেওয়া যায় এমন লক্ষ্য তাপমাত্রা, °C

def cooling_ode(t, T):
    return -K * (T - T_ENV)

print(f"T0 = {T0}°C, T_env = {T_ENV}°C, k = {K} /মিনিট, T_safe = {T_SAFE}°C")
print(f"গভর্নিং সমীকরণ: dT/dt = -k(T - T_env)")

# ========== ধাপ ২: RK4 ODE সলভার + বাইসেকশন রুট-ফাইন্ডার (L56-এর টুলকিট) ==========
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

H_MAIN = 1.0
N_STEPS_MAIN = 90
main_trace = solve_ode(cooling_ode, 0.0, T0, H_MAIN, N_STEPS_MAIN)

print("RK4 তাপমাত্রা-ট্রেস (প্রতি ১০ মিনিটে):")
for t, T in main_trace:
    if round(t) % 10 == 0:
        print(f"  t={t:5.1f} মিনিট  T={T:7.3f}°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, T0
    for _ in range(n_steps):
        y = rk4_step(cooling_ode, t, y, h_actual)
        t += h_actual
    return y

def bisection_root(g, a, b, tol=1e-8, 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 g_safe(t):
    return solve_to_time(t, H_MAIN) - T_SAFE

t_safe_numeric, iters_root = bisection_root(g_safe, 0.0, 120.0, tol=1e-8)
print(f"\nবাইসেকশন রুট-ফাইন্ড: t_safe (RK4, h={H_MAIN}) = {t_safe_numeric:.6f} মিনিট, {iters_root} ইটারেশনে")
print(f"যাচাই: T(t_safe) = {solve_to_time(t_safe_numeric, H_MAIN):.6f}°C")

# ========== ধাপ ৩: ক্লোজড-ফর্ম ভ্যালিডেশন ==========
import math

def T_exact(t):
    return T_ENV + (T0 - T_ENV) * math.exp(-K * t)

def t_exact_safe():
    return -math.log((T_SAFE - T_ENV) / (T0 - T_ENV)) / K

t_safe_exact = t_exact_safe()
error_time = abs(t_safe_numeric - t_safe_exact)

print(f"ক্লোজড-ফর্ম সঠিক t_safe = {t_safe_exact:.6f} মিনিট")
print(f"RK4+বাইসেকশন t_safe     = {t_safe_numeric:.6f} মিনিট")
print(f"সময়ের এবসোলিউট এরর      = {error_time:.8f} মিনিট ({error_time*60:.6f} সেকেন্ড)")

T_final_numeric = main_trace[-1][1]
T_final_exact = T_exact(main_trace[-1][0])
error_temp = abs(T_final_numeric - T_final_exact)
print(f"\nt={main_trace[-1][0]:.1f} মিনিটে: RK4 T={T_final_numeric:.6f}°C, সঠিক T={T_final_exact:.6f}°C, এরর={error_temp:.8f}°C")

# ========== ধাপ ৪: কনভারজেন্স-অর্ডার যাচাই (স্টেপ-সাইজ অর্ধেক করে) ==========
T_TEST_TIME = 60.0
step_sizes = [4.0, 2.0, 1.0, 0.5]

print(f"t={T_TEST_TIME} মিনিটে সঠিক T = {T_exact(T_TEST_TIME):.8f}°C\n")

errors = []
for h in step_sizes:
    n = round(T_TEST_TIME / h)
    trace_h = solve_ode(cooling_ode, 0.0, T0, h, n)
    T_h = trace_h[-1][1]
    err = abs(T_h - T_exact(T_TEST_TIME))
    errors.append(err)
    print(f"  h={h:4.2f}  n_steps={n:4d}  T={T_h:.8f}  এরর={err:.10f}")

print("\nh অর্ধেক করলে এরর-রেশিও (RK4 চতুর্থ-অর্ডারের জন্য প্রত্যাশিত ~16x):")
for i in range(1, len(errors)):
    ratio = errors[i-1] / errors[i]
    print(f"  h={step_sizes[i-1]} -> h={step_sizes[i]}: রেশিও = {ratio:.3f}")

# ========== ধাপ ৫: মন্টে কার্লো সেনসিটিভিটি অ্যানালাইসিস (k-এর অনিশ্চয়তা) ==========
import random
random.seed(7)

K_SD = 0.006   # k-এর অনিশ্চয়তা (স্ট্যান্ডার্ড ডেভিয়েশন), ইলাস্ট্রেটিভ
N_TRIALS = 2000

def t_safe_for_k(k_val):
    def ode(t, T):
        return -k_val * (T - T_ENV)
    def solve_to_time_k(t_target, h):
        n_steps = max(1, round(t_target / h))
        h_actual = t_target / n_steps
        t, y = 0.0, T0
        for _ in range(n_steps):
            y = rk4_step(ode, t, y, h_actual)
            t += h_actual
        return y
    def g(t):
        return solve_to_time_k(t, 1.0) - T_SAFE
    root, _ = bisection_root(g, 0.0, 150.0, tol=1e-6)
    return root

samples = []
for _ in range(N_TRIALS):
    k_sample = max(0.01, random.gauss(K, K_SD))
    samples.append(t_safe_for_k(k_sample))

mean_t = sum(samples) / len(samples)
var_t = sum((x - mean_t)**2 for x in samples) / (len(samples) - 1)
sd_t = math.sqrt(var_t)
samples_sorted = sorted(samples)
p05 = samples_sorted[int(0.05 * len(samples_sorted))]
p95 = samples_sorted[int(0.95 * len(samples_sorted))]

print(f"মন্টে কার্লো সেনসিটিভিটি (n={N_TRIALS}, k ~ N({K}, {K_SD}^2)):")
print(f"  গড় t_safe   = {mean_t:.4f} মিনিট")
print(f"  স্ট্যান্ডার্ড ডেভিয়েশন = {sd_t:.4f} মিনিট")
print(f"  ৫ম পার্সেন্টাইল  = {p05:.4f} মিনিট")
print(f"  ৯৫তম পার্সেন্টাইল = {p95:.4f} মিনিট")
print(f"  নমিনাল (k={K}) t_safe = {t_safe_numeric:.4f} মিনিট")

# ========== ধাপ ৬: চূড়ান্ত সারাংশ ও সুপারিশ -- ধাপ ১-৫-এর প্রকৃত ফলাফল থেকে ==========
print("=" * 60)
print("চূড়ান্ত মূল্যায়ন সারাংশ -- ফার্নেস-যন্ত্রাংশ শীতলীকরণ (হাইপোথেটিক্যাল)")
print("=" * 60)
print(f"১. নিউমেরিক্যাল সমাধান (RK4+বাইসেকশন, h={H_MAIN}): t_safe = {t_safe_numeric:.4f} মিনিট")
print(f"২. ক্লোজড-ফর্ম ভ্যালিডেশন: t_safe = {t_safe_exact:.4f} মিনিট, এরর = {error_time:.8f} মিনিট")
print(f"৩. কনভারজেন্স-অর্ডার: h অর্ধেক করলে এরর ~{errors[0]/errors[1]:.1f}x, ~{errors[1]/errors[2]:.1f}x, ~{errors[2]/errors[3]:.1f}x কমে (তাত্ত্বিক 16x-এর কাছাকাছি)")
print(f"৪. সেনসিটিভিটি (k-এর অনিশ্চয়তা): গড় = {mean_t:.3f} মিনিট, sd = {sd_t:.3f} মিনিট, [P5,P95] = [{p05:.3f}, {p95:.3f}]")

recommended = mean_t + 2 * sd_t
print(f"\nরক্ষণশীল অপেক্ষা-সময়ের সুপারিশ (গড় + ২×sd, ~৯৭.৫% নিরাপত্তা মার্জিন): {recommended:.2f} মিনিট")

if error_time < 0.01 and 16.0 * 0.5 <= (errors[1]/errors[2]) <= 16.0 * 1.5:
    accuracy_verdict = "নিউমেরিক্যাল মেথড অত্যন্ত সঠিক ও তাত্ত্বিক অর্ডার অনুযায়ী আচরণ করছে"
else:
    accuracy_verdict = "নিউমেরিক্যাল মেথডের অ্যাকুরেসি বা কনভারজেন্স-অর্ডার পুনরায় যাচাই প্রয়োজন"

if sd_t / mean_t > 0.10:
    reliability_verdict = "কিন্তু k-এর অনিশ্চয়তার কারণে t_safe-এর প্রকৃত মান যথেষ্ট অনিশ্চিত (relative sd > 10%)"
else:
    reliability_verdict = "এবং k-এর অনিশ্চয়তা t_safe-এর উপর সীমিত প্রভাব ফেলে (relative sd <= 10%)"

print(f"\nমেথড-অ্যাকুরেসি সিদ্ধান্ত: {accuracy_verdict}")
print(f"নির্ভরযোগ্যতা সিদ্ধান্ত: {reliability_verdict}")
print(f"\nচূড়ান্ত প্রকৌশলগত সুপারিশ: পরবর্তী অপারেশন t = {recommended:.1f} মিনিটের আগে শিডিউল করবেন না")