微分方程的数值解法—四阶龙格库塔(Runge-Kutta)的推导以及C代码实现

 四阶龙格库塔(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);
}