c++实现 龙格库塔经典4阶算法

1. 龙格-库塔(Runge-Kutta)方法简介

经典四阶法

在各种龙格－库塔法当中有一个方法十分常用，以至于经常被称为“RK4”或者就是“龙格－库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。令初值问题表述如下。

对于该问题的RK4由如下方程给出：

其中

这样，下一个值（yn+1）由现在的值（yn）加上时间间隔（h）和一个估算的斜率的乘积所决定。该斜率是以下斜率的加权平均：

k1是时间段开始时的斜率；

k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；

k3也是中点的斜率，但是这次采用斜率k2决定y值；

k4是时间段终点的斜率，其y值用k3决定。

当四个斜率取平均时，中点的斜率有更大的权值：

RK4法是四阶方法，也就是说每步的误差是h阶，而总积累误差为h阶。

注意上述公式对于标量或者向量函数（y可以是向量）都适用。

以上内容来自链接:https://www.jianshu.com/p/e4aa9a688959

---

2. 原文章中使用的是matlab直接转换出来的c语言代码,拷贝到编译器后会出现很多error,所以索性用c++重写,顺便做个记录

1. 代码参照原文,对其中的一些变量和算法中的过程做了简单注释
2. 对容易崩溃的指针操作点进行了使用前检查

#include <iostream>

double dEPS = 1.0e-16;	//精度
//---------------------网上的例子--------------------------
double* Func_aa(double* R, int R_size, double* const dRdt) {
	std::fill(dRdt, dRdt + R_size, 0.0);
	if (R_size <= 0) return dRdt;
	dRdt[0] = R[1];      /*y0'=y1*/
	dRdt[1] = -R[0];     /*y1'=y0*/
	dRdt[2] = -R[2];     /*y2'=y2*/
	return dRdt;
}


//---------------------下面是龙格库塔函数---------------------
using Func_xx = double* (*)(double* R, int R_size, double* const dRdt);

/* 参数说明
* param_num:微分方程中未知数的个数
* step:		积分的步数(包括这一步)
* start_t	积分的起点
* steplength:积分的步长
* R0:		方程Func_xx的初始值入参
* z[]:		二维数组,体积为 param_num * step , 返回step个积分点上的param_num个函数值
*/
void RKT(const int& param_num, int step, double start_t, double steplength, 
double* R0, double* z, Func_xx myFunc)
{
	double* inparam_buf = nullptr;		//输入参数的缓存
	double* inparam = nullptr;			//当前输入参数
	double* valueOut = nullptr;		//myFunc的结果

	if (param_num < 1)
	{
		std::cout << "error:积分参数个数为0" << std::endl;
		return;
	}

	int bufsize = param_num * sizeof(double);
	valueOut = (double*)malloc(bufsize);
	inparam_buf = (double*)malloc(bufsize);
	inparam = (double*)malloc(bufsize);

	if (inparam_buf == nullptr|| inparam == nullptr|| valueOut == nullptr)
	{
		std::cout << "内存分配失败" << std::endl;
		return;
	}

	double a[3] = { 0 };
	// 初始化龙格库塔公式中的系数	
	a[0] = steplength / 2.0;		//这是用来求k2的时候用的 1/2
	a[1] = steplength / 2.0;		//k2的系数 
	a[2] = steplength;				//k3的系数	
	a[3] = steplength;				//k4的系数

	memcpy(inparam, R0, bufsize);	// 输入初始化
	// 二维数组初始化
	for(size_t i = 0 ; i < param_num ; ++i)
	{
		z[i] = inparam[i];	//使用第一组参数初始化第一次的结果
	}
	for (size_t w = 1; w < step; ++w)
	{
		myFunc(inparam, param_num, valueOut);
		for (size_t p_index=0; p_index < param_num; ++p_index )
		{
			inparam_buf[p_index] = inparam[p_index];	//初值
		}
		for (size_t f = 0; f < 3; ++f)
		{
			for (size_t yi = 0 ; yi < param_num ; ++yi)
			{
				//  y(tn) = yn; //yn + +a[i] * k1
				inparam[yi] = z[(w - 1) * param_num + yi] + a[f] * valueOut[yi];  
				//yn+1 = yn + h/6(k1+k2+k3+k4)
				inparam_buf[yi] = inparam_buf[yi] + a[f + 1] * valueOut[yi] / 3.0;
			}
			myFunc(inparam, param_num, valueOut);
		}

		for (size_t i = 0 ; i < param_num ; ++i)
		{
			//最后的加上k4/6,因为前三次只有k1+k2+k3的合
			inparam[i] = inparam_buf[i]  + steplength * valueOut[i] / 6.0;
		}
		for (size_t i= 0 ; i<param_num ; ++i)
		{
			z[w * param_num + i] = inparam[i];			//保存每次的结果
		}
		start_t = start_t + steplength;					// tn = tn-1 + steplength
	}
	free(inparam_buf);
	free(inparam);
	free(valueOut);
	valueOut = NULL;
	inparam_buf = NULL;
	inparam = NULL;
}

// main.cpp
int main()
{
    double end_t = 0.0001;
    double start_t = 0.0;
    const int step = 512;
    const int paramnum = 3;
    
    double R[3] = { 2,3,4 };
    double dRdt[3] = { 0 };
    double* z = new double[step * paramnum];
    RKT(paramnum, step, start_t, end_t / step, R, z, Func_aa);
    double zz[step][paramnum] = { 0 };
    memcpy(zz, z, sizeof(double) * paramnum*step); 
	//zz[i]表示每次的结果数组一共512组数据
	
    system("pause");
}

结果如图: