Επίλυση διαφορικής εξίσωσης με τη μέθοδο Runge-Kutta στη C++

Η Runge-Kutta είναι μία προσεγγιστική μέθοδος επίλυσης διαφορικών εξισώσεων. Στην ουσία πρόκειται για μια μεγάλη οικογένεια μεθόδων με πολλές παραλλαγές.

Η μέθοδος είναι περισσότερο ακριβής από τη μέθοδο Euler (Δείτε Υπολογιστικές μέθοδοι για επίλυση διαφορικών εξισώσεων) και στην πράξη χρησιμοποιείται πολύ περισσότερο.

Εδώ θα εξετάσουμε ένα απλό παράδειγμα της κλασικής μεθόδου Runge-Kutta (RK4). Έστω η διαφορική εξίσωση:

$$ \frac{\mathrm{d}\,y}{\mathrm{d}\,x}= x\,y $$

με αρχικές συνθήκες x=0, y=1. Να υπολογιστεί πίνακας τιμών της συνάρτησης y=f(x) στο διάστημα [0,1].

Η μέθοδος στηρίζεται στον αναδρομικό τύπο:

$$ \begin{align} y_{n+1} &= y_n + \tfrac{1}{6} h\left(k_1 + 2k_2 + 2k_3 + k_4 \right)\\ x_{n+1} &= x_n + h \\ \end{align} $$

όπου:

$$ \begin{align} k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + \tfrac{1}{2}h , y_n + \tfrac{h}{2} k_1) \\ k_3 &= f(x_n + \tfrac{1}{2}h , y_n + \tfrac{h}{2} k_2) \\ k_4 &= f(x_n + h , y_n + hk_3). \end{align} $$

Ένα απλό πρόγραμμα που επιλύει τη διαφορική εξίσωση είναι:

#include<iostream>
#include<stdio.h>
#include<math.h>

#define NMAX 500

using namespace std;

double f1(double x, double y)
{
    // return dy/dx
    return x*y;
}

int main()
{
    int i,n;
    double a, b, h;
    double x[NMAX],  y[NMAX];
    double k1, k2, k3, k4;

    /* limits */
    cout << "Please enter lower limit, a= " ;
    cin >> a;
    cout << "Please enter lower limit, b= " ;
    cin >> b;

    /* discretization */
    cout << "Please enter number of steps: (<=500!) ";
    cin >> n;
    if(n > NMAX)
    {
        cout << "Please try again with fewer steps\n";
        return 0;
    }
    /* initial conditions */
    cout << "Please enter initial value y[0]= ";
    cin >> y[0];
  
    h  = (b - a) / n;      /* step */
    x[0] = a;  /* initial value of x */

    /* loop calculations */
    for(i=0; i < n; i++)
    {
        k1 = h * f1(x[i],y[i]);
        k2 = h * f1(x[i] + h/2.0, y[i] + 0.5 * k1);
        k3 = h * f1(x[i] + h/2.0, y[i] + 0.5 * k2);
        k4 = h * f1(x[i] + h, y[i] + k3);
        y[i+1] = y[i] + (k1 + 2.0*k2 +2.0*k3 + k4 ) / 6.0;
        x[i+1] = a + h*(i+1);
    }

    /* print resulted table */
    for (i=0; i <= n; i++)
       printf("%5d %10.6f %10.6f\n", i, x[i], y[i]);
  
    return 0;
}

Συνδεθείτε για περισσότερες δυνατότητες αλληλεπίδρασης,
σχολιασμοί, εξωτερικοί σύνδεσμοι, βοήθεια, ψηφοφορίες, αρχεία, κτλ.

Εκπαιδευτικό υλικό από τον Αθανάσιο Σταυρακούδη σας παρέχετε κάτω από την άδεια Creative Commons Attribution-NonCommercial-ShareAlike 4.0 License.
Σας παρακαλώ να ενημερωθείτε για κάποιους επιπλέον περιορισμούς
http://stavrakoudis.econ.uoi.gr/stavrakoudis/?iid=401.

Athanassios Stavrakoudis, Sat 23 Nov 2024