গসিয়ান কোয়াড্রেচার

১ · প্রশ্ন — পয়েন্ট বাছাই কি স্বাধীন করা যায়?

ট্র্যাপিজয়ডাল ও সিম্পসনস রুল (L23-L24) — এবং নিউটন-কোটস পরিবারের সব মেথড — একটি সাধারণ বৈশিষ্ট্য শেয়ার করে: এগুলো সবসময় সমান-দূরত্বের পয়েন্টে ফাংশনের মান ব্যবহার করে, শুধু সেই পয়েন্টগুলোর ওয়েট (গুণক) পরিবর্তন করে নির্ভুলতা বাড়ায়। গসিয়ান কোয়াড্রেচার একটি ভিন্ন প্রশ্ন করে — যদি পয়েন্টের অবস্থানও স্বাধীনভাবে বেছে নেওয়ার স্বাধীনতা থাকে (শুধু ওয়েট নয়), তাহলে একটি নির্দিষ্ট সংখ্যক পয়েন্টে সর্বোচ্চ কত নির্ভুলতা অর্জন করা সম্ভব?

উত্তরটি চমকপ্রদ — nটি পয়েন্ট ব্যবহার করে সমান-দূরত্বের নিউটন-কোটস মেথড সাধারণত ডিগ্রি n-1 (বা প্রতিসাম্যের সুবিধায় সিম্পসনসের মতো n) পলিনোমিয়াল পর্যন্ত নিখুঁত। কিন্তু যদি পয়েন্টের অবস্থান স্বাধীনভাবে বেছে নেওয়া যায়, একই nটি পয়েন্ট দিয়ে ডিগ্রি 2n-1 পর্যন্ত পলিনোমিয়াল নিখুঁতভাবে ইন্টিগ্রেট করা সম্ভব — প্রায় দ্বিগুণ ক্ষমতা!

২ · গস-লিজেন্ড্র নোড ও ওয়েট

[-1, 1] ইন্টারভালের জন্য গস-লিজেন্ড্র মেথডের সর্বোত্তম নোড (xᵢ) ও ওয়েট (wᵢ) লিজেন্ড্র পলিনোমিয়ালের মূল থেকে গাণিতিকভাবে ডেরাইভ করা হয় — এগুলো সুপরিচিত ধ্রুবক, টেক্সটবুকে হার্ডকোড করা মান হিসেবে ব্যবহার করা স্ট্যান্ডার্ড প্র্যাকটিস।

২-পয়েন্ট গস-লিজেন্ড্র (ডিগ্রি-৩ পলিনোমিয়াল পর্যন্ত নিখুঁত):

৩-পয়েন্ট গস-লিজেন্ড্র (ডিগ্রি-৫ পলিনোমিয়াল পর্যন্ত নিখুঁত):

যেকোনো নির্বিচার [a, b] ইন্টারভালে প্রয়োগ করতে একটি সরল লিনিয়ার ম্যাপিং ব্যবহার করা হয়:

৩ · সত্যিকারের ডেমো — অ্যাকুরেসি-পার-ইভ্যালুয়েশন তুলনা

L23-L24-এর একই পরীক্ষার ফাংশন (∫₀^π sin(x) dx = 2) ব্যবহার করে, এবার গস-লিজেন্ড্র মেথড হাতে-লেখা কোড দিয়ে গণনা করে সিম্পসনস রুলের সাথে সমান সংখ্যক ফাংশন-ইভ্যালুয়েশনে তুলনা করা হলো:

import math

def f(x):
    return math.sin(x)

a, b = 0.0, math.pi
true_val = 2.0

def gauss2(f, a, b):
    # স্ট্যান্ডার্ড 2-পয়েন্ট গস-লিজেন্ড্র নোড ও ওয়েট, [-1,1]-এর জন্য
    nodes   = [-1/math.sqrt(3), 1/math.sqrt(3)]
    weights = [1.0, 1.0]
    mid, half = (a+b)/2, (b-a)/2
    total = sum(wi * f(mid + half*xi) for xi, wi in zip(nodes, weights))
    return half * total

def gauss3(f, a, b):
    # স্ট্যান্ডার্ড 3-পয়েন্ট গস-লিজেন্ড্র নোড ও ওয়েট, [-1,1]-এর জন্য
    nodes   = [-math.sqrt(3/5), 0.0, math.sqrt(3/5)]
    weights = [5/9, 8/9, 5/9]
    mid, half = (a+b)/2, (b-a)/2
    total = sum(wi * f(mid + half*xi) for xi, wi in zip(nodes, weights))
    return half * total

def simpson(f, a, b, n):
    h = (b-a)/n
    total = f(a) + f(b)
    for i in range(1, n):
        x = a + i*h
        total += (4 if i % 2 != 0 else 2) * f(x)
    return total * h / 3

g2, g3 = gauss2(f, a, b), gauss3(f, a, b)
err_g2, err_g3 = abs(g2-true_val), abs(g3-true_val)
print(f"2-point Gauss-Legendre (2 evals): approx={g2:.10f}  error={err_g2:.3e}")
print(f"3-point Gauss-Legendre (3 evals): approx={g3:.10f}  error={err_g3:.3e}")
print()

s2, s4 = simpson(f, a, b, 2), simpson(f, a, b, 4)
err_s2, err_s4 = abs(s2-true_val), abs(s4-true_val)
print(f"Simpson n=2 (3 evals): approx={s2:.10f}  error={err_s2:.3e}")
print(f"Simpson n=4 (5 evals): approx={s4:.10f}  error={err_s4:.3e}")