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++重写,顺便做个记录

代码参照原文,对其中的一些变量和算法中的过程做了简单注释

对容易崩溃的指针操作点进行了使用前检查

#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/2a[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] * k1inparam[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");
}

结果如图:

c++实现龙格库塔经典4阶算法相关推荐

[计算机数值分析]四阶龙格-库塔经典格式解常微分方程的初值问题
龙格-库塔方法的设计思想: 四阶龙格-库塔方法的经典格式: 程序设计框图: 例:设取步长h=0.2,从x=0到x=1用四阶经典格式解决以下常微分方程的初值问题. 运行示例: 源码: #include& ...
4阶经典龙格库塔格式
%用途:4阶经典龙格库塔格式解常微分方程组y'=f(x,y),y(x0)=y0 %格式:[x,y]=marunge4s(dyfun,xspan,y0,h) %dyfun为向量函数f(x,y),xspa ...
Ode45以及龙格-库塔算法
Ode45及龙格-库塔算法 ODE45 ODE45函数描述: 求解非刚性微分方程 - 中阶方 ODE是Matlab专门用于解微分方程的功能函数.该求解器有变步长(variable-step)和定步长( ...
捷联惯导基础知识之一（姿态更新的毕卡算法、龙格库塔算法、及精确数值解法）
一.介绍旋转矢量以及由旋转矢量更新四元数等效旋转矢量微分方程(Bortz方程): 经过理论和近似推导: 单位四元数与旋转矢量关系: 单位四元数与三角函数的关系: 四元数更新方程:利用旋转矢量,考虑了 ...
四阶龙格库塔c语言,四阶龙格库塔算法的C语言实现
解微分方程 2001年3月焦作大学学报 JOURNALOFJIAOZUOUNIVERSITY№ 1Mar．2001第1期四阶龙格一库塔算法的C语言实现毋玉芝 (焦作财会学校) 摘要本文叙述了四阶龙 ...
四阶龙格库塔算法用MATLAB写
四阶龙格-库塔算法可以使用 MATLAB 进行编写.您可以使用 MATLAB 的数值解法工具箱来解决常微分方程组,并使用相应的函数(例如 ode45)来实现四阶龙格-库塔算法.在编写代码时,您需要根据 ...
matlab:使用4阶龙格库塔方法求解常微分方程组
%书籍:常用数值算法及其matlab实现 %第10章常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
matlab：使用4阶龙格库塔方法求微分方程组的值
参考书目:<数值分析(matlab版)> 作者:周璐翻译 %调用龙格库塔方法求解微分方程组,P399,9.7 微分方程组 function [t,z] = my_rk4b(f, t0, ...
龙格-库塔(Runge-Kutta)方法数学原理及实现
参考:https://blog.csdn.net/u013007900/article/details/45922331 龙格-库塔(Runge-Kutta)方法龙格-库塔(Runge-Kutta) ...
龙格-库塔（Runge-Kutta）法解微分方程
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法.由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂.该算法是构建在数学支持的基础之上的. 对于一阶精度的欧 ...

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

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

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

经典四阶法

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

c++实现龙格库塔经典4阶算法相关推荐

最新文章

热门文章

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

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

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

经典四阶法

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

c++实现 龙格库塔经典4阶算法相关推荐

最新文章

热门文章

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

c++实现龙格库塔经典4阶算法相关推荐