四阶Runge-Kutta数学原理
给定常微分方程初值问题
{ dy dx = f ( x , y ) , a ≤ x ≤ b y ( a ) = α h = b − a N \left\{ \begin{matrix} \frac{\text{dy}}{\text{dx}} = f(x,y),a \leq x \leq b \\ \ y(a) = \alpha\text{\;\;\quad\;h} = \frac{b - a}{N} \\ \end{matrix} \right.\ {dxdy=f(x,y),a≤x≤b y(a)=αh=Nb−a
记 x n = a + n ⋅ h x_{n} = a + n \cdot h xn=a+n⋅h, n = 0 , 1 , ⋯ , N n = 0,1,\cdots,N n=0,1,⋯,N,利用四阶龙格—库塔方法
K 1 = h f ( x n , y n ) K_{1} = hf\left( x_{n},y_{n} \right) K1=hf(xn,yn)
K 2 = h f ( x n + h 2 , y n + K 1 2 ) K_{2} = hf\left( x_{n} + \frac{h}{2},y_{n} + \frac{K_{1}}{2} \right) K2=hf(xn+2h,yn+2K1)
K 3 = h f ( x n + h 2 , y n + K 2 2 ) K_{3} = hf\left( x_{n} + \frac{h}{2},y_{n} + \frac{K_{2}}{2} \right) K3=hf(xn+2h,yn+2K2)
K 4 = h f ( x n + h , y n + K 3 ) K_{4} = hf\left( x_{n} + h,y_{n} + K_{3} \right) K4=hf(xn+h,yn+K3)
y n + 1 = y n + 1 6 ( K 1 + 2 K 2 + 2 K 3 + K 4 ) y_{n + 1} = y_{n} + \frac{1}{6}(K_{1} + 2K_{2} + 2K_{3} + K_{4}) yn+1=yn+61(K1+2K2+2K3+K4)
n = 0 , 1 , ⋯ , N − 1 n = 0,1,\cdots,N - 1 n=0,1,⋯,N−1
可逐次求出微分方程初值问题的数值解 y n y_{n} yn, n = 1 , 2 , ⋯ , N n = 1,2,\cdots,N n=1,2,⋯,N。
程序流程
核心代码
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
int n;
double a, b, fa;
double f(double x, double y) { return -y * y; }
double f_(double x) { return 1.0 / (x + 1.0); }
int main() {
scanf("%lf%lf%lf%d", &a, &b, &fa, &n);
double x = a, y = fa, h = (b - a) / n;
for (int i = 1; i <= n; i++) {
double k1 = h * f(x, y);
double k2 = h * f(x + h / 2, y + k1 / 2);
double k3 = h * f(x + h / 2, y + k2 / 2);
double k4 = h * f(x + h, y + k3);
x += h;
y += 1.0 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
printf("%.2lf\t%lf\t%.2lf\n", x, y, fabs(f_(x) - y));
}
return 0;
}