#59. Całkowanie numeryczne – metoda trapezów

Dzisiejszy temat nie należy do najłatwiejszych. Całkowanie numeryczne i metoda trapezów brzmi groźnie, ale postaram się przekazać Ci jak najłatwiejszym językiem informacje na temat tego, w jaki sposób całkować metodą trapezów.

Ale zanim zaczniemy mam dla Ciebie ciekawą informację. Czy wiesz, co łączy młodego chłopaka, który dorabia jako kelner, nie rozumie matematyki na studiach i wszystko, co zapisane jest w matematycznych książkach, uważa za bełkot i niemal siwieje, ucząc się matematyki z całkami i rachunkiem różniczkowym?

Na pierwszy rzut oka wydawać by się mogło, że nic. A co jeśli powiem Ci, że tym młodym facetem był Izaac Newton, który jest uważany za niezależnego twórcę rachunku różniczkowego!

Nadal uważasz, że matematyka nie jest dla Ciebie? Wstydź się!

Nie jest to też pierwszy wpis na temat całkowania na blogu. Popełniłem jeden, podobny wpis na temat całkowania metodą prostokątów.

Czas zacząć!

Metoda trapezów — wprowadzenie

Wspomniana metoda prostokątów nie należy do najlepszych sposób na całkowanie. Dużo bardziej wydajna jest metoda trapezów. Nie będę omawiał tu podstaw całkowania, ponieważ zrobiłem to we wspomnianym wyżej artykule. O ile go znasz, kojarzysz zapewne rysunek, w którym pokazany był wykres podzielony na wąskie prostokąty. To kluczowy element całkowania, ponieważ używamy go do tego, by obliczyć pole powierzchni pod wykresem funkcji.

Przyglądając się wykresom poniżej, możemy zauważyć, że w przypadku metody prostokątów nie do końca dokładnie odwzorowana jest powierzchnia pod krzywą. Im większych używamy prostokątów, tym większy błąd w obliczeniach. Z trapezami jest nieco inaczej, ponieważ bardziej dopasowują się do krzywizny.

Wykresy te są mocno uproszczone, ale widać jak na tacy, że metoda trapezów jest dokładniejsza. 

Każdy trapez ma wysokość równą dx, a podstawy równe są długości w krańcach podziałów. W przypadku powyżej posiadamy przedział całkowania od a do b, który podzielony został na kilka trapezów lub prostokątów. Naszym zadaniem jest stworzenie takiego algorytmu, który na wejście otrzyma:

• funkcję
• przedział całkowania — wartość początkowa i końcowa
• ilość trapezów, na jaki ma zostać podzielony przedział całkowania

Na wyjściu otrzymamy wartość całki w przybliżeniu.

Wzór

Czas na trochę więcej matematyki. Definicja całki oznaczonej Riemana mówi nam, że wartość całki równa jest sumie pól obszarów pod wykresem krzywej w zadanym przedziale całkowania. Sumę możemy uzyskać tak, jak wspominałem, dzieląc obszar na mniejsze podobszary.

Zobaczmy na poniższe wzory, które bardziej dotykają szkoły podstawowej niż studiów. Po pierwsze szerokość trapezu, którą obliczamy, będzie później naszą wysokością we wzorze na pole trapezu. Z racji podzielenia obszaru na równej szerokości trapezy możemy zastosować wzór, gdzie odejmujemy prawą część obszaru (b) od lewej (a) i dzielimy przez ilość trapezów.

Następnie wysokości trapezów, które później we wzorze numer trzy będą podstawami. Wzór jest intuicyjny, aby obliczyć bok lewy i prawy trapezu i-tego w kolejności, musimy od boku po lewej stronie odjąć numer trapezu przemnożony przez jego szerokość wyliczaną w pierwszym wzorze.

Na końcu pozostaje nam obliczyć pole trapezu ze wzoru, jaki każdy z nas powinien znać. Sumy podstaw dzielone przez dwa i na końcu pomnożone przez wysokość. 

Po takich obliczeniach pozostaje nam zsumować wszystkie otrzymane pola. Opis jest dość długi, ale jak rozrysujesz sobie wszystko na kartce papieru, to okaże się, że jest to mega proste!

Metoda trapezów — Implementacja

Dobrze, teorię mamy za sobą. Nie wprowadzałem celowo skomplikowanych wzorów, bo nie o to chodzi. Programowanie ma być przyjemnością, nawet gdy dotykamy „potwornej” matematyki. Praktyka pokaże nam, że całki nie są takie strasznie złe.

Rozwiązanie stworzę w pięciu językach. W języku Python 3 przekazujemy o jeden parametr więcej na start, ponieważ zastosujemy dla odmiany funkcję lambda.

Zdefiniujemy metodę o nazwie trapezoidal_rule, która przyjmie trzy lub cztery parametry — funkcję, początek i koniec zakresu oraz liczbę trapezów, na które podzielimy naszą figurę pod wykresem.

Na początku metody deklarujemy zmienną integral, która przechowa wartość całki. W kolejnym kroku obliczymy szerokość trapezu (dx). Czas na najważniejszą część funkcji całkującej. Pętla, w której to pętli korzystamy ze wzoru numer dwa, gdzie obliczamy wartości podstaw górnej i dolnej.

Na końcu sumujemy każde pole trapezu do zmiennej integral, używając operatora +=. Pozostaje nam zwrócenie wyniku, a następnie wywołanie funkcji z parametrami i cieszenie się wynikiem.

Dla przykładu funkcji x + 2 wynikiem będzie całka o wartości przybliżonej równej 2.5. 

Python 3

def trapezoidal_rule(fun, a, b, n):
  integral = 0

  dx  = (b-a)/n

  for i in range(n):
    fa = a + dx * i
    fb = a + dx * (i + 1)

    integral += (fun(fa) + fun(fb)) / 2 * dx
  return integral

integral = trapezoidal_rule (lambda x: x+2, 0.0, 1.0, 100)
print(integral)

Java

package com.company;

public class Main {

    private static double fun(double x) {
        return x+2;
    }

    public static double trapezoidal_rule(double a, double b, double n ) {
        double integral = 0;

        double dx = (b - a)/n;

        for (int i = 0; i < n; i++)
        {
            double fa = a + dx * i;
            double fb = a + dx * (i + 1);

            integral += (fun(fa) + fun(fb)) / 2 * dx;
        }
        return integral;
    }

    public static void main(String[] args) {
        double integral = trapezoidal_rule(0, 1, 100);
        System.out.println(Double.toString(integral));
    }
}

C#

using System;

public class Program
{
	private static double fun(double x) 
	{
        return x+2;
    }
	
	public static void Main()
	{
		double integral = trapezoidal_rule(0, 1, 100);
  		Console.WriteLine(integral.ToString());
	}
	
	
    public static double trapezoidal_rule(double a, double b, double n ) 
	{
        double integral = 0;

        double dx = (b - a)/n;

        for (int i = 0; i < n; i++)
        {
            double fa = a + dx * i;
            double fb = a + dx * (i + 1);

            integral += (fun(fa) + fun(fb)) / 2 * dx;
        }
        return integral;
    }
}

C++

#include <iostream>
using namespace std;
 
double fun(double x) 
{
   return x+2;
}

double trapezoidal_rule(double a, double b, double n ) 
{
    double integral = 0;

    double dx = (b - a)/n;

    for (int i = 0; i < n; i++)
    {
        double fa = a + dx * i;
        double fb = a + dx * (i + 1);

        integral += (fun(fa) + fun(fb)) / 2 * dx;
    }
        
    return integral;
}


int main() {
  double integral = trapezoidal_rule(0, 1, 100);
  cout << integral;
  return 0;
}
   

JavaScript

var fun = function(x) {         
    return x + 2;
};

var trapezoidal_rule = function(a, b, n) {         
    var integral = 0;

    var dx = (b - a)/n;

    for (var i = 0; i < n; i++)
    {
       var fa = a + dx * i;
       var fb = a + dx * (i + 1);
       integral += (fun(fa) + fun(fb)) / 2 * dx;
    }
   return integral;
};

Metoda trapezów — podsumowanie

Analiza matematyczna i świat całek to ogromne zagadnienia.  Z jednej strony straszne (przez pryzmat studiów technicznych), z drugiej jednak fascynujący.

Ja Tobie pokazałem tylko niewielki fragment, który dotyka całek oznaczonych. Mimo wszystko metoda trapezów ukazuje piękno i wielkie możliwości, jakie stoją za rachunkiem różniczkowym. Sam nie uważam się za perfekcyjnie władającego tą częścią matematyki i nawet powiem więcej — uważam, że znam analizę matematyczną tylko w niewielkim stopniu. Może czas nadrobić?

60 wpis na blogu to rozmowa na temat bootcampów!

Newsletter

Nie przegap i dołącz już dziś do 838 osób będących w tym Newsletter! Otrzymuj co niedzielę o godzinie 20 listę kilku ciekawych tematów, które miałem okazję obserwować w mijającym tygodniu.

Tematy będą głównie techniczne, ale czasami pojawi się coś, co może wprowadzi Cię w stan zadumy i zmusi do dyskusji w szerszym gronie. Zero spamu!