Całkowanie numeryczne - metoda interpolacji wielomianem drugiego stopnia

Autor podstrony: Krzysztof Zajączkowski

Stronę tą wyświetlono już: 3179 razy

Całkowanie numeryczne metodą interpolacji drugiego stopnia umożliwia dokładniejsze przeprowadzenie obliczeń przybliżonej całki z zadanej funkcji f(x) w przedziale od a do b. Do obliczenia całki z takiej interpolowanej funkcji potrzebny będzie następujący wzór:

[1]

Zapis wyrażenia w formacie TeX-a:

\int_{a}^{b}(a_w\cdot x^2+b_w\cdot x+c_w)\,dx=\left[\frac{1}{3}\cdot a_w\cdot x^3+\frac{1}{2}\cdot b_w\cdot x^2+c_w\cdot x\right]_{a}^{b}=\frac{1}{3}\cdot a_w\cdot b^3+\frac{1}{2}\cdot b_w\cdot b^2+c_w\cdot b - \frac{1}{3}\cdot a_w\cdot a^3-\frac{1}{2}\cdot b_w\cdot a^2-c_w\cdot a

Do interpolacji liniowej konieczne będzie wykorzystanie trzech punktów, które złożą się na następujący układ równań:

[2]

Zapis wyrażenia w formacie TeX-a:

\begin{cases} a_w\cdot {x_1}^2 +b_w\cdot x_1+c_w = y_1 \\ a_w\cdot {x_2}^2 +b_w\cdot x_2+c_w = y_2 \\ a_w\cdot {x_3}^2 +b_w\cdot x_3+c_w = y_3 \end{cases}

Do wyznaczenia współczynników a_w, b_w i c_w wielomianu interpolującego wykorzystane zostaną wzory Cramera, które opisane zostały na stronie Matematyka → Równania liniowe → Rozwiązywanie układów równań za pomocą wzorów Cramera oraz wzór na obliczenie wyznacznika macierzy 3×3 z strony Matematyka → Macierze → Wyznacznik macierzy.

W językach programowania takich jak C++ do wykorzystania tego algorytmu konieczne jest zaimplementowanie algorytmu ONP, który powstał jako modyfikacja algorytmu notacji polskiej Jana Łukasiewicza. W przypadku języków skryptowych takich jak PHP czy Python można posłużyć się funkcją eval, której użycie ze względów bezpieczeństwa jest niekiedy niewskazane.

Oto przykład prostego programu obliczającego tą metodą całkę podanej na wejście funkcji napisany w Pythonie:

#!/usr/bin/env python from math import * class Point2D: def __init__(self, x, y): self.x = x self.y = y def __str__(self): return "x = {x}; y = {y}".format(x = self.x, y = self.y) class Matrix3x3: def __init__(self, point1, point2, point3): self.point1 = point1 self.point2 = point2 self.point3 = point3 def getDet(self): return self.point1.x * self.point1.x * self.point2.x + self.point2.x * self.point2.x * self.point3.x + self.point3.x * self.point3.x * self.point1.x - self.point2.x * self.point3.x * self.point3.x -self.point3.x * self.point1.x * self.point1.x - self.point1.x * self.point2.x * self.point2.x def getDetA(self): return self.point1.y * self.point2.x + self.point2.y * self.point3.x + self.point3.y * self.point1.x - self.point2.x * self.point3.y - self.point3.x * self.point1.y - self.point1.x * self.point2.y def getDetB(self): return self.point1.x * self.point1.x * self.point2.y + self.point2.x * self.point2.x * self.point3.y + self.point3.x * self.point3.x * self.point1.y - self.point2.y * self.point3.x * self.point3.x - self.point3.y * self.point1.x * self.point1.x - self.point1.y * self.point2.x * self.point2.x def getDetC(self): return self.point1.x * self.point1.x * self.point2.x * self.point3.y + self.point2.x * self.point2.x * self.point3.x * self.point1.y + self.point3.x * self.point3.x * self.point1.x * self.point2.y - self.point1.y * self.point2.x * self.point3.x * self.point3.x - self.point2.y * self.point3.x * self.point1.x * self.point1.x - self.point3.y * self.point1.x * self.point2.x * self.point2.x; def getABC(self): w = self.getDet() return [self.getDetA() / w, self.getDetB() / w, self.getDetC() / w] def integrate(function, a, b, i): dx = (b - a) / i integr = 0 for x in range(i): x = x * dx + a xa = x xb = x + dx p1 = Point2D(x, eval(function)) x += dx / 2 p2 = Point2D(x, eval(function)) x += dx / 2 p3 = Point2D(x, eval(function)) interpolation = Matrix3x3(p1, p2, p3) abc = interpolation.getABC() integr += abc[0] / 3 * (xb * xb * xb - xa * xa * xa) + abc[1] / 2 * (xb * xb - xa * xa) + abc[2] * (xb - xa) return integr def main(args): function = input("Podaj funkcję: ") a = float(input("Podaj a: ")) b = float(input("Podaj b: ")) i = int(input("Podaj i: ")) print("Całka z funkcji {function} w przedziale od {a} do {b} jest równa {integrate}".format(function = function, a = a, b = b, integrate = integrate(function, a, b, i))) return 0 if __name__ == '__main__': import sys sys.exit(main(sys.argv))

Przykład działania dla f(x) = x², a = 0, b = 1, oraz i = 10:

Podaj funkcję: x**2
Podaj a: 0
Podaj b: 1
Podaj i: 10
Całka z funkcji x**2 w przedziale od 0.0 do 1.0 jest równa 0.333333333333517
Aby kontynuować, naciśnij dowolny klawisz . . .

Jak widać metoda ta sprawdza się bardzo dobrze, zwłaszcza gdy obliczana funkcja jest wielomianem drugiego stopnia, warto jednak sprawdzić, jak zadziała np. dla funkcji f(x)=sin(x) w przedziale od 0 do 3.14:

Podaj funkcję: sin(x)
Podaj a: 0
Podaj b: 3.14
Podaj i: 10
Całka z funkcji sin(x) w przedziale od 0.0 do 3.14 jest równa 2.000005502397631
Aby kontynuować, naciśnij dowolny klawisz . . .

Uzyskany powyżej wynik jest zawyżony, gdyż rzeczywista wartość tej całki powinna wynosić 1.999998731727541, ale różnica jest rzędu 6.770670088585007·10^-6. Przy 100 przedziałach wynik jest jeszcze dokładniejszy:

Podaj funkcję: sin(x)
Podaj a: 0
Podaj b: 3.14
Podaj i: 100
Całka z funkcji sin(x) w przedziale od 0.0 do 3.14 jest równa 1.9999987323929425
Aby kontynuować, naciśnij dowolny klawisz . . .

Jak widać w tym przypadku różnica wynosi 0 (przynajmniej na tym poziomie zakresu dokładności).