হিট ইকুয়েশনের জন্য ফাইনাইট ডিফারেন্স

import math

alpha = 1.0
Nx = 10                 # স্থানিক ইন্টারভাল সংখ্যা -> ১১টি গ্রিড-বিন্দু
dx = 1.0 / Nx
xs = [i*dx for i in range(Nx+1)]

def initial_condition(x):
    return math.sin(math.pi*x)

def run_ftcs(r, dt, steps):
    u = [initial_condition(x) for x in xs]
    u[0] = 0.0; u[-1] = 0.0
    history = [u[:]]                 # list of lists: প্রতিটি সারি এক সময়-ধাপের প্রোফাইল
    for n in range(steps):
        new = u[:]
        for i in range(1, Nx):
            new[i] = u[i] + r*(u[i+1] - 2*u[i] + u[i-1])
        new[0] = 0.0; new[-1] = 0.0
        u = new
        history.append(u[:])
    return history

r_stable = 0.4                       # CFL-সদৃশ শর্ত: r <= 0.5 (স্থিতিশীল)
dt_stable = r_stable * dx*dx / alpha
print(f"dx={dx}  dt={dt_stable:.6f}  r={r_stable}")

steps = 5
hist = run_ftcs(r_stable, dt_stable, steps)

print()
print("তাপমাত্রা প্রোফাইলের বিবর্তন (স্থিতিশীল, r=0.4):")
print("ধাপ  " + "".join(f"x={x:.1f}".rjust(9) for x in xs))
for n, row in enumerate(hist):
    print(f"{n:4d} " + "".join(f"{v:9.5f}" for v in row))

t_final = steps*dt_stable
print()
print(f"t_final = {t_final:.4f}")
print(f"{'x':>6} {'FD u(x,t)':>12} {'নির্ভুল':>12} {'|এরর|':>10}")
for i, x in enumerate(xs):
    exact = math.sin(math.pi*x) * math.exp(-math.pi**2 * alpha * t_final)
    fd = hist[-1][i]
    print(f"{x:6.2f} {fd:12.6f} {exact:12.6f} {abs(fd-exact):10.2e}")

max_err = max(abs(hist[-1][i] - math.sin(math.pi*xs[i])*math.exp(-math.pi**2*alpha*t_final)) for i in range(Nx+1))
print()
print("সর্বোচ্চ |এরর| t_final-এ:", max_err)

def run_ftcs_spike(r, steps):
    u = [0.0]*(Nx+1)
    u[Nx//2] = 1.0     # কেন্দ্রে স্পাইক -- সব ফ্রিকোয়েন্সি মোড উত্তেজিত করে
    history = [u[:]]
    for n in range(steps):
        new = u[:]
        for i in range(1, Nx):
            new[i] = u[i] + r*(u[i+1] - 2*u[i] + u[i-1])
        new[0] = 0.0; new[-1] = 0.0
        u = new
        history.append(u[:])
    return history

print(f"{'r':>6} {'কেন্দ্র বিন্দুর মান, ধাপ 0..8':<70} {'স্থিতিশীল?':>12}")
for r in [0.4, 0.6, 0.9]:
    hist = run_ftcs_spike(r, 8)
    center_vals = [row[Nx//2] for row in hist]
    stable = "হ্যাঁ (r<=0.5)" if r <= 0.5 else "না (r>0.5)"
    vals_str = ", ".join(f"{v:.4f}" for v in center_vals)
    print(f"{r:6.1f} {vals_str:<70} {stable:>12}")