MN_lab_04.py

# Numerical integration - Całkowanie numeryczne
import math
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import pylab

# ======================================================

def f1(x: float) -> float:
    return math.sin(x)


def f2(x: float) -> float:
    return pow(x, 2)+2*x+5


def f3(x: float) -> float:
    return math.exp(x)

# ======================================================

def rectangularIntegration(func, a: float, b: float, n: int):
    s = (b-a)/n
    SUM = 0
    for i in range(0, n):
        SUM += func((a+(i*s))+(1/2)*s)  # f((a+(i*s))+(1/2)*s)
    SUM *= s
    return SUM


def trapezoidalIntegration(func, a: float, b: float, n: int):
    s = (b-a)/n
    SUM = 0
    for i in range(0, n):
        S1 = (((i+1)*s)-(i*s))
        S2 = (func(a+(i*s)))+func(a+((i+1)*s))
        SUM += (S1/2*S2)
    return SUM


def parabolicIntegration(func, a: float, b: float, n:int):
    S1  = 0; S2 = 0;
    s = (b - a) / n;
    for i in range(1,n+1):
        x = a + i * s
        S2 += func(x - s / 2)
        if i < n:
            S1 += func(x)
    SUM = s / 6 * (func(a) + func(b) + 2 * S1 + 4 * S2)
    return SUM

# ======================================================

def GaussLegendr(f, a:float,b:float,n:int,tab):
    SUM=0
    for i in range(0,n):
        t = ((a+b)/2)+(((b-a)/2)*tab['x'][i])
        SUM+=f(t)*tab['A'][i]
    return SUM*((b-a)/2)

        
if __name__ == '__main__':
      tab2={"x":[-math.sqrt(3)/3,math.sqrt(3)/3],"A":[1,1]}
      tab3={"x":[-math.sqrt(3/5),0,math.sqrt(3/5)],"A":[5/9,8/9,5/9]}
      tab4={"x":[-math.sqrt(525+70*math.sqrt(30))/35,-math.sqrt(525-70*math.sqrt(30))/35,math.sqrt(525-70*math.sqrt(30))/35,math.sqrt(525+70*math.sqrt(30))/35],"A":[(18-math.sqrt(30))/36,(18+math.sqrt(30))/36,(18+math.sqrt(30))/36,(18-math.sqrt(30))/36]}
      #print(GaussLegendr(f1,0.5, 2.5,len(tab2['x']),tab2))
      #print(GaussLegendr(f1,0.5, 2.5,len(tab3['x']),tab3))
      #print(GaussLegendr(f1,0.5, 2.5,len(tab4['x']),tab4))

      print("f(x) = sin(x)")
      print("Gauss n=2: ",GaussLegendr(f1,0.5, 2.5,len(tab2['x']),tab2))
      print("Gauss n=3: ",GaussLegendr(f1,0.5, 2.5,len(tab3['x']),tab3))
      print("Gauss n=4: ",GaussLegendr(f1,0.5, 2.5,len(tab4['x']),tab4))
      print("Rect. n=20:", rectangularIntegration(f1, 0.5, 2.5, 20))
      print("Trap. n=20:", trapezoidalIntegration(f1, 0.5, 2.5, 20))
      print("Para. n=20:", parabolicIntegration(f1, 0.5, 2.5, 20))
      print("Wolfram:   ", "1.67873\n\n")
      
      print("f(x) = x^2 + 2x + 5")
      print("Gauss n=2: ",GaussLegendr(f2,0.5, 5,len(tab2['x']),tab2))
      print("Gauss n=3: ",GaussLegendr(f2,0.5, 5,len(tab3['x']),tab3))
      print("Gauss n=4: ",GaussLegendr(f2,0.5, 5,len(tab4['x']),tab4))
      print("Rect. n=20:", rectangularIntegration(f2, 0.5, 5, 20))
      print("Trap. n=20:", trapezoidalIntegration(f2, 0.5, 5, 20))
      print("Para. n=20:", parabolicIntegration(f2, 0.5, 5, 20))
      print("Wolfram:   ", "88.875\n\n")

      print("f(x) = exp(x)")
      print("Gauss n=2: ",GaussLegendr(f3,0.5, 5,len(tab2['x']),tab2))
      print("Gauss n=3: ",GaussLegendr(f3,0.5, 5,len(tab3['x']),tab3))
      print("Gauss n=4: ",GaussLegendr(f3,0.5, 5,len(tab4['x']),tab4))
      print("Rect. n=20:", rectangularIntegration(f3, 0.5, 5, 20))
      print("Trap. n=20:", trapezoidalIntegration(f3, 0.5, 5, 20))
      print("Para. n=20:", parabolicIntegration(f3, 0.5, 5, 20))
      print("Wolfram:   ", "146.764\n\n")