四阶龙格库塔(Runge-Kutta)的推导
对于一阶常微分方程
的解y = y(x),根据微分中值定理,存在一点使得
步长用h表示,即h=x(i+1)-x(i),进一步有
引入平均斜率K
带入y(x(i+1))可表示为
下面以二阶龙格库塔为例,在区间上取两点,以该两点的斜率值K1、K2的加权平均来求取平均斜率K的近似值
K1为Xi点处的切线斜率值
K2为X(i+1)点处的切线斜率值
带入,y(x(i+1))又可表示为
接下来将K1、K2在Xi处泰勒展开
带入上式中,y(x(i+1))进一步可表示为
接下来再将y(x(i+1))在Xi处泰勒展开
最后将两个式子对应项的系数进行对比得
最后另a2 = b21 = 1、c1 = c2 = 0.5,便得到二阶的龙格库塔公式
同理可得四阶经典的龙格库塔公式,二阶是取两点的斜率K1、K2;四阶取4点的斜率K1,K2,K3,K4,推导同二阶的一样。注意每一阶的龙格库塔公式不唯一,c1,c2,a2,b21的取值不同会得到不一样的公式,下面的只是比较经典的,常用的。
C代码实现
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "utf_SLR.h"
#define N 40 //定义点数
typedef void(*FUNTYPE)(double, double*,double *);
void csp_scmove(double t, double* s,double *ds); //微分方程
void ode45(FUNTYPE csp_scmove, double sta, double end, double **s0, double* s1)
{
int i, j;
double h;
double* k1, * k2, * k3, * k4;
double* res1, * res2, * res3, * res4;
double t, * s2, * s3, * s4;
s2 = (double*)malloc(6 * sizeof(double));
s3 = (double*)malloc(6 * sizeof(double));
s4 = (double*)malloc(6 * sizeof(double));
k1 = (double*)malloc(6 * sizeof(double));
k2 = (double*)malloc(6 * sizeof(double));
k3 = (double*)malloc(6 * sizeof(double));
k4 = (double*)malloc(6 * sizeof(double));
res1 = (double*)malloc(6 * sizeof(double));
res2 = (double*)malloc(6 * sizeof(double));
res3 = (double*)malloc(6 * sizeof(double));
res4 = (double*)malloc(6 * sizeof(double));
//}
h = (end - sta) / N;
if (s1 == NULL) {
printf("内存分配不成功!\n");
}
else {
for (i = 0; i < 6; i++) {
s1[i] = s0[i][0];
}
}
for (i = 0; i < N; i++)
{
t = i * h; //x为一系列的点
csp_scmove(t, s1, res1);
for (j = 0; j < 6; j++) {
k1[j] = res1[j];
s2[j] = s1[j] + (h / 2) * k1[j];
}
csp_scmove(t + h / 2, s2, res2);
for (j = 0; j < 6; j++) {
k2[j] = res2[j];
s3[j] = s1[j] + (h / 2) * k2[j];
}
csp_scmove(t + h / 2, s3, res3);
for (j = 0; j < 6; j++) {
k3[j] = res3[j];
s4[j] = s1[j] + h * k3[j];
}
csp_scmove(t + h, s4, res4);
for (j = 0; j < 6; j++) {
k4[j] = res4[j];
}
for (j = 0; j < 6; j++) {
s1[j] = s1[j] + h * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6;
}
}
free(s2); free(s3); free(s4); free(k1); free(k2); free(k3); free(k4);
free(res1); free(res2); free(res3); free(res4);
}