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 0039
//---------------------------

#include <iostream>
#include <windows.h>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy
// dane do zadania
//------------------

// Wysokość, na której
// znajduje się lufa
// działa [m]
double hh = 10;
// Przyspieszenie grawitacyjne
// Ziemi [m/s^2]
double gg = 9.81;
// Prędkość pocisku
// w osi x [m/s]
double Vx = 1000;

// 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

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

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

int main()
{
  // Zmienne

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

  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(4)
       << fixed;

  cout << "Równanie nieliniowe, "
          "rozwiązanie metodą "
          "Steffensena\n"
          "---------------------"
          "-------------------"
          "-----------\n\n";

  // 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 wartość
    // funkcji 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!\n\n";

  cout << "Zasięg strzału    d    = "
       << setw(10) << x
       << " [m]" << endl << endl
       << "Zasięg maksymalny dmax = "
       << setw(10) << dmax
       << " [m]" << endl << endl;
  system("pause");
  return 0;
}

Równanie nieliniowe, rozwiązanie metodą Steffensena
---------------------------------------------------

Zasięg strzału    d    =   819.1981 [m]

Zasięg maksymalny dmax =  1427.8431 [m]

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

# Równanie nieliniowe:
# metoda Steffensena
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0039
#---------------------------

import math

# Tutaj definiujemy
# dane do zadania
#------------------

# Wysokość, na której
# znajduje się lufa
# działa [m]
hh = 10
# Przyspieszenie grawitacyjne
# Ziemi [m/s^2]
gg = 9.81
# Prędkość pocisku
# w osi x [m/s]
Vx = 1000

# Tutaj definiujemy funkcję,
# której pierwiastek jest
# wyliczany
#---------------------------
def fg(x):
    return (hh - gg * x * x / 2 /
           Vx / Vx - math.log(x))

# Tutaj definiujemy
# parametry początkowe

# Dokładność wyznaczania
# pierwiastka
epsx = 1e-7
# Dokładność wyznaczania
# zera funkcji
epsy = 1e-7
# Maksymalna liczba obiegów
n = 64

# Program główny
#---------------
result = False

print("Równanie nieliniowe, "
      "rozwiązanie metodą "
      "Steffensena\n"
      "---------------------"
      "-------------------"
      "-----------\n")

# Obliczamy dmax
dmax = Vx * math.sqrt(2 * hh / gg)

# Punkt startowy
x = dmax

while n > 0:
    n -= 1

    # Obliczamy wartość
    # funkcji w punkcie x
    h = fg(x)

    # Sprawdzamy, czy wartość
    # funkcji jest dostatecznie
    # bliska zeru
    if abs(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 abs(x1 - x) < epsx:
        result = True
        break

    # Kontynuujemy obliczenia

if not result:
    print("Zakończono z błędem!\n")

print(f"Zasięg strzału    d    = "
      f"{x:10.4f} [m]\n\n"
      f"Zasięg maksymalny dmax = "
      f"{dmax:10.4f} [m]\n")
input("Naciśnij Enter...")