求解常微分方程组初值问题的龙格库塔法分析及其C代码.docVIP

求解常微分方程组初值问题的 龙格库塔法分析及其C代码 1、概 述 由高等数学的知识可知，一些特殊类型的常微分方程（组）能够求出给定初始值的解析解，而在科学与工程问题中遇到的常微分方程（组）往往是极其复杂的，要想求得其给定初始值的解析解就变得极其困难，甚至是得不到解析解。尽管如此，在研究实际问题时，往往只需要获得若干点上的近似值就行了，这就说明研究常微分方程（组）的数值解法是很有必要的。求解常微分方程（组）的数值解法有多种，比如欧拉法、龙格库塔法、线性多步法等等，其中龙格库塔法是这几种方法中比较常用的，因为其计算精度较高且具有自启动的特点。 对于用龙格库塔法求解单个常微分方程和求解常微分方程组的思路基本相似（注意一点一个微分方程组是常微分方程组即表明微分方程中的各阶导数都是对同一个变量求导，例如可以把各个量对时间求导得到一个常微分方程组，如果一个微分方程组中的有对不同变量的导数那么这个方程组就成了偏微分方程组），都是根据泰勒展开得到其迭代计算形式，其基本思想都是按照一定原则求取当前点附近一些点的斜率，通过这些斜率的线性组合作为当前点处的斜率，进行递推求解。 2、数学模型 2.1 常微分方程初值问题的数学模型 （1） 为的已知函数，为给定的初始值。 常微分方程的数值解法的任务就是要求出函数值在微分自变量取如下序列： 时的值的近似值，一般情况下都采取等距结点的方式，即：，其中为相邻两结点的距离，称为步长。 2.2 常微分方程组的数学模型 （2） 为的已知函数，为初始值，利用数值解法求解常微分方程组的任务就是求出函数值在微分自变量取如下序列： 时其值的近似值 3 龙格库塔法的迭代形式 龙格库塔法的迭代形式推导在相关的数值分析书籍中均有详细的讲解，这里就不进行推导直接给出常用的四阶龙格库塔法的形式作为例子。 3.1 单个常微分方程的迭代形式 对应于2.1节的数学模型有： 通过上面的迭代形式可知实际上就是如下四个点 （点1） （点2） （点3） （点4） 处关于的斜率。这四个点的确定又是一环扣一环的，比如要求得点2的位置，就必须先求得在点1处的斜率 对应于2.2节的数学模型有： （4） （方程的个数）个，即是要把的斜率一次性算完，才能继续往下递推运算。 在求解常微分方程组初值问题的情况下，可将常微分方程组写成向量的形式，这样更加便于理解，下面用图形来展示整个求解过程： 图1 龙格库塔法的求解 沿着箭头的方向依次解出各个向量，即可依次递推下去 ，下面的表示为点处的斜率，处情形类似。 4 算 例 4.1 问题描述 利用四阶龙格库塔法求解初值问题 取步长为。 4.2 问题分析 为了和第3节的数学模型对应，我们令，，则上面的问题就变为和前面的模型一样的形式: 4.3 C程序代码 #include stdio.h #include math.h #define M 2 // 维数,方程组中方程的个数 // 算斜率的右端函数 void RightSlopeFun( double *pRightK, // 算出来的右端斜率项[out] double iX, // 计算点的分自变量值 double *iY // 计算点的初始值 ) { pRightK[0] = iX*iY[0]-iY[1]; pRightK[1] = (iX+iY[0])/iY[1]; } // 龙格库塔法解算程序 double LgktSoulution( double *pOY, // 返回函数值[out] double iX, // 积分自变量输入值 double *iY, // 初始值 double h, // 单步步长 int Count=1 // 解算步数，默认形参为1，调用时不给此数赋值，只解算一步 ) { double K1[M]; // 第一个斜率 double K2[M]; // 第二个斜率 double K3[M]; // 第三个斜率 double K4[M]; // 第四个斜率 double tempY[M]; // 保存Y的中间值的临时空间 int i; // 计数器 double tempX; // 临时保存x的空间 tempX=iX;