BVP-এর জন্য ফাইনাইট ডিফারেন্স মেথড

import math

# একই BVP: y'' = -y on [0,1],  y(0)=0,  y(1)=2 ; exact y(x) = C*sin(x)
C_exact = 2.0 / math.sin(1.0)

N = 9            # অভ্যন্তরীণ গ্রিড-বিন্দু সংখ্যা
h = 1.0 / (N + 1)
y0, yN1 = 0.0, 2.0   # বাউন্ডারি মান

# সমীকরণ: y_{i-1} + (h^2 - 2)*y_i + y_{i+1} = 0   (i = 1..N)
a = [0.0]*(N+1); b = [0.0]*(N+1); c = [0.0]*(N+1); d = [0.0]*(N+1)
for i in range(1, N+1):
    a[i], b[i], c[i], d[i] = 1.0, -(2.0 - h*h), 1.0, 0.0

# বাউন্ডারি কন্ডিশন ডান-পাশে (rhs) ভাঁজ করা
d[1] -= a[1]*y0
d[N] -= c[N]*yN1

def thomas_solve(a, b, c, d, n):
    cp = [0.0]*(n+1); dp = [0.0]*(n+1)
    cp[1] = c[1]/b[1]; dp[1] = d[1]/b[1]
    for i in range(2, n+1):
        m = b[i] - a[i]*cp[i-1]
        cp[i] = c[i]/m if i < n else 0.0
        dp[i] = (d[i] - a[i]*dp[i-1]) / m
    x = [0.0]*(n+1)
    x[n] = dp[n]
    for i in range(n-1, 0, -1):
        x[i] = dp[i] - cp[i]*x[i+1]
    return x

y = thomas_solve(a, b, c, d, N)
xs = [i*h for i in range(N+2)]
ys = [y0] + [y[i] for i in range(1, N+1)] + [yN1]

print(f"{'x':>6} {'FD y(x)':>14} {'নির্ভুল C*sin(x)':>18} {'|এরর|':>12}")
for i in range(N+2):
    exact = C_exact * math.sin(xs[i])
    print(f"{xs[i]:6.2f} {ys[i]:14.8f} {exact:18.8f} {abs(ys[i]-exact):12.2e}")

# রেসিডুয়াল যাচাই: মূল (আন-ফোল্ডেড) সমীকরণে ফিরে বসিয়ে দেখা
max_res = max(abs(ys[i-1] + b[i]*ys[i] + ys[i+1]) for i in range(1, N+1))
print()
print("সর্বোচ্চ |রেসিডুয়াল| (হওয়া উচিত ~0):", max_res)

# কনভারজেন্স-অর্ডার যাচাই: N দ্বিগুণ করলে এরর কীভাবে কমে (x=0.5-এ)
print()
print("গ্রিড সূক্ষ্ম করলে x=0.5-এ এরর (O(h^2) কনভারজেন্স যাচাই):")
for Ntest in [9, 19, 39]:
    h2 = 1.0/(Ntest+1)
    a2=[0.0]*(Ntest+1); b2=[0.0]*(Ntest+1); c2=[0.0]*(Ntest+1); d2=[0.0]*(Ntest+1)
    for i in range(1, Ntest+1):
        a2[i], b2[i], c2[i], d2[i] = 1.0, -(2.0-h2*h2), 1.0, 0.0
    d2[Ntest] -= c2[Ntest]*yN1
    y2 = thomas_solve(a2, b2, c2, d2, Ntest)
    idx = round(0.5/h2)
    y_at_half = y2[idx]
    err = abs(y_at_half - C_exact*math.sin(idx*h2))
    print(f"N={Ntest:3d}  h={h2:.5f}  y(0.5)≈{y_at_half:.8f}  এরর={err:.3e}")