Równania

W równaniu nieliniowym zmienna x jest argumentem funkcji nieliniowej. Ogólnie możemy zapisać równanie w sposób następujący:

Po przekształceniu otrzymujemy:

Tworzymy nową funkcję:

I mamy równanie:

Problem sprowadza się zatem do znalezienia pierwiastków funkcji h(x). Każdy pierwiastek tej funkcji jest jednocześnie rozwiązaniem równania wyjściowego. Pierwiastki możemy szukać jedną z opisanych wcześniej metod pod warunkiem, że znamy przedział, w którym znajduje się pierwiastek. I tutaj mamy podstawową trudność. W przypadku ogólnym nie istnieje algorytm, który znajduje wszystkie pierwiastki równania metodami przybliżonymi. Możemy jedynie szukać pierwiastków w pewnym przedziale argumentów funkcji. Nie jest to zadanie proste. Zwykle jednak przy rozwiązywaniu problemów technicznych możemy określić przedział, w którym interesuje nas znalezienie rozwiązania. W takim przypadku możemy wykorzystać jedną z metod obliczania przybliżonego pierwiastka funkcji, np. metodę Steffensena, ponieważ wymaga ona tylko jednego punktu startowego.

Dla przykładu rozwiążmy następujące zadanie z fizyki.

Na wysokości h = 10 m znajduje się działo, którego lufa umieszczona jest poziomo. Z działa zostaje wystrzelony pocisk, który opuszcza lufę z prędkością V = 1000 m/s. Przed działem znajduje się teren, który wznosi się wg funkcji:

Na pocisk działa tylko siła ciążenia ziemskiego. Pominąć tarcie powietrza. Należy wyznaczyć odległość d od końca lufy działa, w jakiej upadnie pocisk.

Zobaczmy na rysunek poglądowy:

Najpierw rozważmy lot pocisku tak, jakby teren był płaski. Gdy kula opuści lufę armaty, to jej prędkość rozkłada się na dwie składowe:

Składowa pozioma V_x jest stała, ponieważ na pocisk nie działa opór powietrza, który pomijamy w zadaniu w celu uproszczenia.

Składowa pionowa V_y pochodzi od przyciągania ziemskiego i wyraża się wzorem:

gdzie: g  – ziemskie przyspieszenie grawitacyjne równe 9,81 m/s2 t – czas trwania ruchu

Jeśli teren jest płaski, to ruch pocisku będzie trwał tak długo, jak spada on w pionie z wysokości h, zatem:

Z kolei w tym czasie pocisk przesunie się w poziomie na odległość d_max:

Jeśli teren się wznosi, to pocisk zderzy się z nim wcześniej w odległości d, która znajduje się w przedziale od 0 do d_max. Musimy teraz wyznaczyć funkcję wysokości pocisku w zależności od x. Oznaczmy tę funkcję jako f(x). W poziomie położenie x pocisk osiąga po czasie t_x:

Z kolei po czasie t_x pocisk znajduje się na wysokości h_x:

Łączymy obie zależności i otrzymujemy funkcję wysokości pocisku f(x):

Wysokość terenu określa funkcja g(x):

Aby pocisk uderzył w teren, musi się znaleźć na wysokości terenu, zatem:

Tworzymy różnicę funkcji:

Szukamy pierwiastka funkcji fg(x), który jest równy odległości d:

Przedział poszukiwań to (0,d_max]. Za punkt startowy wybierzemy środek tego przedziału, czyli d_max / 2.

Poniższy program pokazuje w sposób przykładowy, jak znaleźć odległość d.

// Równanie nieliniowe - metoda Steffensena
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//------------------------------------------

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy dane do zadania
//----------------------------------
double hh = 10;   // Wysokość, na której znajduje się lufa działa
double gg = 9.81; // Przyspieszenie grawitacyjne Ziemi
double Vx = 1000; // Prędkość pocisku w osi x

// Tutaj definiujemy funkcję, której pierwiastek jest wyliczany
//-------------------------------------------------------------

double fg(double x)
{
  return hh - gg * x * x / 2 / Vx / Vx - log(x) ;
}

// Tutaj definiujemy parametry początkowe

double epsx = 1e-7; // Dokładność wyznaczania pierwiastka
double epsy = 1e-7; // Dokładność wyznaczania zera
int n       = 64;   // Maksymalna liczba obiegów

// Program główny
//---------------

int main()
{
    // Zmienne

    double x,x1,h,g,dmax;
    bool result = false;

    setlocale(LC_ALL,"");

    cout << fixed;
    cout << "Równanie nieliniowe, rozwiązanie metodą Steffensena" << endl
         << "---------------------------------------------------" << endl << endl;

    // Obliczamy dmax

    dmax = Vx * sqrt(2 * hh / gg);

    // Punkt startowy

    x = dmax;

    while(--n)
    {
       // Obliczamy wartość funkcji w punkcie x

       h = fg(x);

       // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

       if(fabs(h) < epsy)
       {
           result = true;
           break;
       }

       // Obliczamy różnicę dzieloną

       g = (fg(x+h) - h) / h;

       // Zapamiętujemy bieżące przybliżenie

       x1 = x;

       // Obliczamy kolejne przybliżenie

       x -= h/g;

       // Sprawdzamy, czy odległość pomiędzy dwoma ostatnimi przybliżeniami
       // jest mniejsza od założonej dokładności

       if(fabs(x1 - x) < epsx)
       {
           result = true;
           break;
       }

       // Kontynuujemy obliczenia
    }

    if(!result) cout << "Zakończono z błędem!" << endl << endl;

    cout << "Zasięg strzału    d    = " << setw(18) << x << endl << endl;
    cout << "Zasięg maksymalny dmax = " << setw(18) << dmax << endl << endl;
    return 0;
}

Poeksperymentuj z tym programem. Zmieniaj dane wejściowe, funkcję terenu, punkt startowy, dokładności epsilon, itp.