PROGRAMOWANIE – PRZYKŁADY
IDE C++ WinBGI C# Wpadki

Przykład C++

Liniowa aproksymacja średniokwadratowa Algorytm Moduł obliczeniowy Rzeczywisty i ekranowy układ kartezjański Zakres współrzędnych Program w C++ Poprzedni przykład Następny przykład Program w Visual C# Kontakt

Liniowa aproksymacja średniokwadratowa

Jednym z zagadnień pojawiających się przy rozwiązy­waniu wielu problemów fizy­cznych jest odtwo­rzenie funkcji na podstawie pomiarów jej wartości w skoń­czonej liczbie punktów. Zazwyczaj charakter funkcji f opisu­jącej dane zjawisko nie jest dokładnie znany, a wyniki pomiarów są obarczone błędami. W takich przypad­kach można jedynie zadowolić się rozwią­zaniem przybli­żonym, znajdując funkcję g w ustalonej klasie funkcji, np. w postaci wielo­mianu, która aproksy­muje (przybliża) funkcję f.

W najprostszym przypadku poszukiwaną funkcją aproksy­mującą jest funkcja liniowa (wielo­mian stopnia co najwyżej pierwszego), a metoda jej wyzna­czania jest oparta na minima­lizacji odchy­lenia średniokwa­dratowego. W celu uściślenia problemu zakładamy, że w punktach x1,x2,...,xn dane są wartości y1,y2,...,yn (dokładne lub przybli­żone) funkcji f:

Zagadnienie liniowej aproksymacji średniokwa­dratowej, znane w niektórych środo­wiskach zaintere­sowań pod nazwą regresji liniowej, polega na znale­zieniu funkcji g postaci

takiej, aby suma kwadratów odchyleń

była jak najmniejsza. Odchylenia funkcji aproksymu­jącej od punktów (xk,yk) są mierzone wzdłuż osi y (rys.).

Algorytm

Wyrażenie S traktowane jako funkcja dwóch zmiennych rzeczy­wistych a i b jest wielo­mianem stopnia drugiego o warto­ściach nieuje­mnych, ma więc minimum. Aby wyznaczyć wartości a i b, dla których S osiąga to minimum, przyrównu­jemy do zera pochodne cząstkowe S względem a i b. W rezul­tacie otrzymu­jemy układ dwóch równań liniowych o dwóch niewia­domych a i b:

Rozwiązujemy go w taki sposób, że najpierw obliczamy cztery sumy:

a następnie wartości a i b, korzystając ze wzorów Cramera wyraża­jących rozwiązanie układu poprzez wartości jego odpowie­dnich wyzna­czników:

Oczywiście do tych samych wzorów końcowych można dojść bez korzy­stania z wzorów Cramera, rozwią­zując powyższy układ równań metodą podsta­wiania lub przeci­wnych współczyn­ników. Układ ma dokładnie jedno rozwią­zanie (jest nieoso­bliwy) wtedy i tylko wtedy, gdy mianownik w powyższych wyraże­niach jest różny od zera. Łatwo zauważyć, że metoda wymaga co najmniej dwóch punktów.

W obliczeniach komputerowych operacje arytme­tyczne na liczbach rzeczywi­stych (typu floatdouble w C++) nie są dokładne, toteż należy unikać dzielenia nie tylko przez zero, lecz także przez wartości bliskie zero. Z tego względu w dalszych rozważa­niach przyjmiemy, że wartość bezwzględna miano­wnika w powyż­szych wzorach nie może być niższa od ½.10-8. Użytko­wnik może określić inną minimalną wartość dodatnią jako istotnie różną od zera.

Moduł obliczeniowy

Opisany w języku matematycznym algorytm wyzna­czania prostej, która najlepiej przybliża dany zbiór punktów w sensie aproksy­macji średniokwa­dratowej, wyrazimy w postaci modułu C++. Dostępna w nim funkcja, którą nazwiemy linapr, będzie zawierać sześć argu­mentów (parame­trów): n – liczba punktów, x i y – tablice współrzę­dnych tych punktów, a i b – współczyn­niki szukanej prostej, eps – minimalna wartość dodatnia trakto­wana jako istotnie różna od zera. Naturalnie pierwszy argument funkcji będzie przeka­zywany przez wartość. Dwa kolejne jej argu­menty reprezen­tujące tablice jednowy­miarowe możemy zapisać w czytelnej postaci tabli­cowej double x[]double y[]. W rzeczywi­stości będą nimi przeka­zywane przez wartość wskaźniki na pierwsze elementy tych tablic zapisywane w mniej przej­rzysty, ale równo­ważny sposób jako double *xdouble *y. Argu­menty a i b reprezen­tujące wynik obliczeń najwygo­dniej jest przekazać przez referencję, poprze­dzając je symbolem &. Agument eps określa minimalną wartość bezwzględną miano­wnika, która ma być traktowana jako różna od zera (wartości niższe wskazują na układ osobliwy – zob. wzory powyżej). Jest to argument domyślny. Jeśli nie jest podany, jego wartością jest stała linaprEPS. Plik nagłów­kowy linapr.h rozpatry­wanego modułu może mieć postać:

// linapr.h - liniowa aproksymacja średniokwadratowa
// -----------------------------------------------------------
// Parametry funkcji linapr:
//   n     - liczba punktów
//   x, y  - współrzędne punktów (tablice n-elementowe)
//   a, b  - współczynniki prostej y = ax+b (wynik obliczeń)
//   eps   - minimalna wartość dodatnia istotnie różna od zera
// Wartość zwrotna linapr:
//   true  - pomyślny przebieg obliczeń
//   false - obliczenia przerwane (układ osobliwy)
// -----------------------------------------------------------

#ifndef H_LINAPR
#define H_LINAPR

const double linaprEPS = 0.5e-8;

bool linapr(int n, double x[], double y[], double &a, double &b, double eps = linaprEPS);

#endif // H_LINAPR

Drugi plik modułu o nazwie linapr.cpp zawierający kod źródłowy funkcji linapr zapowie­dzianej w pliku nagłów­kowym może wyglądać następu­jąco:

#include "linapr.h"

bool linapr(int n, double x[], double y[], double &a, double &b, double eps)
{
    double sx = 0, sy = 0, sxx = 0, sxy = 0, det;
    for (int k = 0; k < n; k++)
    {
        sx += x[k];
        sy += y[k];
        sxx += x[k] * x[k];
        sxy += x[k] * y[k];
    }
    if ((det = sxx * n - sx * sx) > -eps && det < eps)
        return false;
    a = (sxy * n - sx * sy) / det;
    b = (sxx * sy - sx * sxy) / det;
    return true;
}

Rzeczywisty i ekranowy układ kartezjański

Zazwyczaj rysunek jest definiowany w rzeczy­wistym układzie współrzę­dnych kartezjań­skich. W celu zobrazo­wania go na ekranie urządzenia grafi­cznego konieczne jest przejście od tego układu rzeczywi­stego do układu ekranowego. Operacja ta, zwana okienkowaniem (ang. windowing), dotyczy dwóch obszarów prosto­kątnych o bokach równole­głych do osi współrzę­dnych: okna w układzie rzeczy­wistym (rys. poniżej z lewej) i widoku (ang. viewport) w układzie ekra­nowym (rys. poniżej z prawej). Często widok nie zajmuje całego ekranu. Na przykład może być otoczony ramką lub może być oprócz niego wyświe­tlana dodatkowa infor­macja tekstowa, tabelka lub inny element graficzny.

Załóżmy, że okno w układzie rzeczywistym jest określone warto­ściami grani­cznymi xRmin, xRmax, yRminyRmax, zaś widok w układzie ekranowym warto­ściami xEmin, xEmax, yEminyEmax. Wówczas odwzoro­wanie opisujące przejście od współrzę­dnych (x,y) okna do współrzę­dnych (xE,yE) widoku ma postać:

gdzie sx i sy są czynnikami skalują­cymi spełnia­jącymi zależności

Zakres współrzędnych

Zobrazowanie punktów na ekranie monitora kompute­rowego wymaga określenia mini­malnych i maksy­malnych wartości współrzę­dnych okna w układzie rzeczy­wistym i widoku w układzie ekranowym. Wyzna­czanie ograniczeń xRmin, xRmax, yRminyRmax okna rozpoczy­namy od iteracji przebiega­jącej wszystkie punkty:

xRmin = xRmax = x[0];
yRmin = yRmax = y[0];
for (int k = 1; k < n; k++)
{
    if (x[k] < xRmin)
        xRmin = x[k];
    else if (x[k] > xRmax)
        xRmax = x[k];
    if (y[k] < yRmin)
        yRmin = y[k];
    else if (y[k] > yRmax)
        yRmax = y[k];
}

Kierując się względami estetycznymi, przesuwamy nieco (np. o 1/20 szerokości okna) lewy brzeg okna w lewo i jego prawy brzeg w prawo, aby końce odcinka prostej aproksy­mującej punkty wystawały poza ich zakres wzdłuż osi x. Korekty może również wymagać dolne i górne ograni­czenie okna, bowiem końce odcinka prostej nie powinny znajdować się poza zakresem okna wzdłuż osi y:

d = 0.05 * (xRmax - xRmin);
xRmin -= d;
xRmax += d;
if ((yPoc = a * xRmin + b) < yRmin)
    yRmin = yPoc;
else if (yPoc > yRmax)
    yRmax = yPoc;
if ((yKon = a * xRmax + b) < yRmin)
    yRmin = yKon;
else if (yKon > yRmax)
    yRmax = yKon;

Nie jest to jednak koniec korekty granic okna, ponieważ nie możemy zaakce­ptować jego zerowej wysokości, która wiodłaby do nadmiaru przy obli­czaniu czynnika skalują­cego sy. Przypadek taki zachodzi wtedy, gdy wszystkie punkty mają tę samą współ­rzędną y. Jeśli jest ona istotnie różna od zera, rozsze­rzamy okno w górę i w dół o jej wartość bezwzględną, w przeci­wnym razie rozsze­rzamy okno w obu tych kierun­kach o stałą linaprEPS wyrażającą najmniejszą wartość dodatnią trakto­waną jako istotnie różną od zera:

if (yRmax - yRmin < linaprEPS)
{
    if (yRmax > linaprEPS)
        d = yRmax;
    else if (yRmax < -linaprEPS)
        d = -yRmax;
    else
        d = linaprEPS;
    yRmax += d;
    yRmin -= d;
}

Ograniczenia xEmin, xEmax, yEminyEmax widoku możemy łatwo określić, znając rozmiar obszaru robo­czego. Mając ponownie na uwadze względy estety­czne i odwrotny kierunek osi y okna grafi­cznego, brzegi widoku możemy wyzna­czyć następująco:

xEmin = getmaxx() / 20;
xEmax = getmaxx() - xEmin;
yEmax = getmaxy() / 20;
yEmin = getmaxy() - yEmax;

Program w C++

Zaprezentowany na poniższym listingu program wczytuje współrzędne punktów z pliku teksto­wego do dwóch tablic dynami­cznych, wyznacza prostą przybliża­jącą te punkty w sensie aproksy­macji średniokwa­dratowej i wyświetla w oknie konsoli jej współczyn­niki. Następnie znajduje okno w układzie rzeczy­wistym zawierające wczytane punkty, a po otwarciu trybu grafi­cznego wyznacza widok w układzie ekranowym i czynniki skalujące oraz rysuje odcinek prostej i aproksy­mowane przez nią punkty, po czym zwalnia pamięć zajmowaną przez obie tablice dynamiczne. Na koniec, gdy użytko­wnik naciśnie dowolny znak na klawia­turze, program zamyka tryb graficzny i kończy działanie.

#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <graphics.h>
#include "linapr.h"

using namespace std;

int n = 0;                          // Liczba punktów
double *x = NULL, *y = NULL;        // Tablice współrzędnych punktów
double a, b;                        // Współczynniki prostej y = ax+b
double xRmin, xRmax, yRmin, yRmax;  // Parametry układu rzeczywistego
int xEmin, xEmax, yEmin, yEmax;     // Parametry układu ekranowego
double yPoc, yKon;                  // Wartości końcowe y = ax+b
double sx, sy;                      // Czynniki skalujące

void error(string err)
{
    if (x != NULL) delete[] x;
    if (y != NULL) delete[] y;
    cout << err << endl;
    getchar();
}

bool czytajdane()
{
    cout << "Nazwa pliku: ";
    string s;
    cin >> s;
    ifstream plik(s.c_str());
    if (!plik)
    {
        error("Nieudane otwarcie pliku.");
        return false;
    }
    if (!(plik >> n))
    {
        error("Nieudany odczyt liczby ca\x88kowitej.");
        return false;
    }
    if (n < 2)
    {
        error("Wymagane co najmniej dwa punkty.");
        return false;
    }
    x = new double[n];
    y = new double[n];
    for (int k = 0; k < n; k++)
        if (!(plik >> x[k] >> y[k]))
        {
            error("Nieudany odczyt liczby rzeczywistej.");
            return false;
        }
    plik.close();
    return true;
}

void min_max()
{
    xRmin = xRmax = x[0];
    yRmin = yRmax = y[0];
    for (int k = 1; k < n; k++)
    {
        if (x[k] < xRmin)
            xRmin = x[k];
        else if (x[k] > xRmax)
            xRmax = x[k];
        if (y[k] < yRmin)
            yRmin = y[k];
        else if (y[k] > yRmax)
            yRmax = y[k];
    }
}

void korekta()
{
    double d = 0.05 * (xRmax - xRmin);
    xRmin -= d;
    xRmax += d;
    if ((yPoc = a * xRmin + b) < yRmin)
        yRmin = yPoc;
    else if (yPoc > yRmax)
        yRmax = yPoc;
    if ((yKon = a * xRmax + b) < yRmin)
        yRmin = yKon;
    else if (yKon > yRmax)
        yRmax = yKon;
    if (yRmax - yRmin < linaprEPS)
    {
        if (yRmax > linaprEPS)
            d = yRmax;
        else if (yRmax < -linaprEPS)
            d = -yRmax;
        else
            d = linaprEPS;
        yRmax += d;
        yRmin -= d;
    }
}

int xE(double x)
{
    return int(sx * (x - xRmin)) + xEmin;
}

int yE(double y)
{
    return int(sy * (y - yRmin)) + yEmin;
}

void rysuj()
{
    xEmin = getmaxx() / 20;
    xEmax = getmaxx() - xEmin;
    yEmax = getmaxy() / 20;
    yEmin = getmaxy() - yEmax;
    sx = (xEmax - xEmin) / (xRmax - xRmin);
    sy = (yEmax - yEmin) / (yRmax - yRmin);
    line(xEmin, yE(yPoc), xEmax, yE(yKon));
    setcolor(LIGHTCYAN);
    for (int k = 0; k < n; k++)
    {
        int xp = xE(x[k]), yp = yE(y[k]);
        circle(xp, yp, 3);
        floodfill(xp, yp, LIGHTCYAN);
    }
}

int main()
{
    if (!czytajdane())
        return 1;
    if (!linapr(n, x, y, a, b))
    {
        error("Nieprawid\x88owe dane (uk\x88" "ad osobliwy).");
        return 2;
    }
    cout << "y = ax + b\n--------------\n";
    cout.setf(ios::fixed, ios::floatfield);
    cout << "a = " << setprecision(5) << a << endl;
    cout << "b = " << setprecision(5) << b << endl;
    min_max();
    korekta();
    initwindow(720, 480, "Aproksymacja średniokwadratowa");
    rysuj();
    delete[] x;
    delete[] y;
    getch();
    closegraph();
    return 0;
}

Program składa się z szeregu funkcji pełniących rozmaite zadania. Funkcja czytajdane wczytuje liczbę punktów i ich współ­rzędne ze strumienia C++ kojarzo­nego z plikiem tekstowym o podanej z klawia­tury nazwie do dwóch tablic dynami­cznych. Plik powinien zawierać, oddzielone spacjami lub znakami końca wiersza, liczbę całkowitą i ciąg liczb rzeczywi­stych, po dwie dla każdego punktu. W przy­padku braku żądanego pliku, nieudanego odczytu danych lub niedozwo­lonej liczby punktów program kończy działanie, wypisując stosowny komunikat i zwracając kod powrotu 1. Wykonanie programu może również zakończyć się kodem powrotu 2 oznacza­jącym, że dane zostały wczytane, lecz obli­czenia realizo­wane przez funkcję linapr zakoń­czyły się fiaskiem z powodu próby dzielenia przez zero (układ równań osobliwy). Wystąpienie błędu powoduje wywołanie funkcji error, która zwalnia pamięć obu tablic dynami­cznych, jeśli uprzednio została im przydzie­lona, wypisuje komunikat informu­jący o zaistniałej sytuacji i wstrzy­muje zamknięcie okna konsoli i zakoń­czenie programu do momentu naciśnięcia klawisza. Zauważmy, że do zatrzy­mania programu posłużyła funkcja getchar, gdyż funkcji getch z biblio­teki WinBGI powinno się używać tylko w trybie graficznym.

W przypadku pomyślnego wczytania danych i wykonania obliczeń program wypisuje współczyn­niki znale­zionej prostej, a następnie wyznacza za pomocą funkcji min_maxkorekta okno w układzie rzeczy­wistym zawierające podane punkty, po czym przechodzi w tryb graficzny oraz rysuje te punkty i aproksy­mującą je prostą, korzystając z funkcji rysuj. Jak łatwo zauważyć, funkcja ta najpierw wyznacza widok w układzie ekranowym i czynniki skalujące używane przez funkcje xEyE przekształ­cające współ­rzędne układu rzeczywi­stego na ekranowy. Dopiero po tych wstępnych oblicze­niach rysuje odcinek prostej i punkty w postaci maleńkich okręgów koloru LIGHTCYAN wypełnionych kolorem WHITE. Na koniec program zwalnia pamięć tablic dynami­cznych, a gdy użytkownik naciśnie dowolny znak na klawia­turze, zamyka tryb graficzny i kończy działanie kodem powrotu 0.

A oto wynik wykonania programu dla przykła­dowego pliku Punkty.txt zawiera­jącego współ­rzędne dziesięciu punktów:


Opracowanie przykładu: styczeń 2019