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

c = 1.0
Nx = 20
dx = 1.0 / Nx
xs = [i*dx for i in range(Nx+1)]

def initial_pulse(x):
    # x=0.3-কেন্দ্রিক ত্রিভুজাকার পালস, অর্ধ-প্রস্থ 0.1, উচ্চতা 1
    if 0.2 <= x <= 0.4:
        return 1.0 - abs(x - 0.3)/0.1
    return 0.0

r = 1.0                       # Courant সংখ্যা: c*dt/dx  (CFL শর্ত: r <= 1)
dt = r*dx/c
print(f"dx={dx}  dt={dt}  Courant r={r}")

def run_wave(r, steps):
    u_prev = [initial_pulse(x) for x in xs]
    u_prev[0] = 0.0; u_prev[-1] = 0.0
    u_curr = [0.0]*(Nx+1)                       # প্রথম ধাপ: শুরুর বেগ 0 ধরে বিশেষ সূত্র
    for i in range(1, Nx):
        u_curr[i] = u_prev[i] + 0.5*(r**2)*(u_prev[i+1] - 2*u_prev[i] + u_prev[i-1])
    u_curr[0] = 0.0; u_curr[-1] = 0.0
    history = [u_prev[:], u_curr[:]]            # list of lists: প্রতি সারি এক সময়-ধাপের প্রোফাইল
    for n in range(steps - 1):
        u_next = [0.0]*(Nx+1)
        for i in range(1, Nx):
            u_next[i] = 2*u_curr[i] - u_prev[i] + (r**2)*(u_curr[i+1] - 2*u_curr[i] + u_curr[i-1])
        u_next[0] = 0.0; u_next[-1] = 0.0
        u_prev, u_curr = u_curr, u_next
        history.append(u_curr[:])
    return history

steps = 8
hist = run_wave(r, steps)

print()
print("তরঙ্গের বিবর্তন u(x,t), ধাপ 0..8 (Courant r=1):")
print("ধাপ " + "".join(f"x={x:.2f}".rjust(7) for x in xs))
for n, row in enumerate(hist):
    print(f"{n:3d} " + "".join(f"{v:7.3f}" for v in row))

print()
print("প্রতি ধাপে সর্বোচ্চ বিন্দুর অবস্থান (peak tracking):")
for n, row in enumerate(hist):
    idx = max(range(len(row)), key=lambda i: row[i])
    print(f"ধাপ {n}: idx={idx}  x={xs[idx]:.2f}  height={row[idx]:.4f}")