Phương pháp Runge-Kutta bậc 4 cho phương trình vi phân

Giới thiệu

Phương pháp Runge-Kutta được dùng để giải phương trình vi phân thông thường (ODE). Nó sử dụng hàm dy/dx cho x và y, và cũng cần giá trị ban đầu của y, tức là y(0). Nó tìm giá trị xấp xỉ của y cho x cho trước. Để giải ODE, ta phải tuân theo các công thức sau:

Ở đây h là chiều cao của khoảng chia.

Lưu ý: Từ các công thức này, ta có thể dùng hai k1 và k2 đầu tiên để tìm nghiệm Runge-Kutta bậc 2 cho ODE.

Đầu vào và Đầu ra

Đầu vào: Giá trị x0 và f(x0): 0 và 0; giá trị x = 0.4; giá trị h = 0.1
Đầu ra: Nghiệm của phương trình vi phân: 0.0213594

Thuật toán

rungeKutta(x0, y0, x, h)
Đầu vào: Giá trị x và y ban đầu, giá trị x đích, và chiều cao khoảng h.
Đầu ra: Giá trị y tại x.

Begin
    iteration := (x – x0)/h
    y = y0
    for i := 1 to iteration, do
        k1 := h*f(x0, y)
        k2 := h*f((x0 + h/2), (y + k1/2))
        k3 := h*f((x0 + h/2), (y + k2/2))
        k4 := h*f((x0 + h), (y + k3))
        y := y + (1/6)*(k1 + 2k2 + 2k3 + k4)
        x0 := x0 + h
    done
    return y
End

Ví dụ (C++)

#include <iostream>
using namespace std;

double diffOfy(double x, double y) {
    return ((x*x)+(y*y)); // hàm x^2 + y^2
}

double rk4thOrder(double x0, double y0, double x, double h) {
    int iteration = int((x - x0)/h);  // tính số lần lặp
    double k1, k2, k3, k4;
    double y = y0;  // ban đầu y là f(x0)
    for(int i = 1; i <= iteration; i++) {
        k1 = h*diffOfy(x0, y);
        k2 = h*diffOfy((x0+h/2), (y+k1/2));
        k3 = h*diffOfy((x0+h/2), (y+k2/2));
        k4 = h*diffOfy((x0+h), (y+k3));
        y += double((1.0/6.0)*(k1+2*k2+2*k3+k4));  // cập nhật y bằng delta y
        x0 += h;  // cập nhật x0 thêm h
    }
    return y;  // giá trị f(x)
}

int main() {
    double x0, y0, x, h;
    cout << "Enter x0 and f(x0): "; cin >> x0 >> y0;
    cout << "Enter x: "; cin >> x;
    cout << "Enter h: "; cin >> h;
    cout << "Answer of differential equation: " << rk4thOrder(x0, y0, x, h);
}

Kết quả

Enter x0 and f(x0): 0 0
Enter x: 0.4
Enter h: 0.1
Answer of differential equation: 0.0213594