欢迎来到天天文库
浏览记录
ID:13159716
大小:372.50 KB
页数:8页
时间:2018-07-21
《求解常微分方程组初值问题的龙格库塔法分析及其c代码》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、求解常微分方程组初值问题的龙格库塔法分析及其C代码1、概述由高等数学的知识可知,一些特殊类型的常微分方程(组)能够求出给定初始值的解析解,而在科学与工程问题中遇到的常微分方程(组)往往是极其复杂的,要想求得其给定初始值的解析解就变得极其困难,甚至是得不到解析解。尽管如此,在研究实际问题时,往往只需要获得若干点上的近似值就行了,这就说明研究常微分方程(组)的数值解法是很有必要的。求解常微分方程(组)的数值解法有多种,比如欧拉法、龙格库塔法、线性多步法等等,其中龙格库塔法是这几种方法中比较常用的,因为其计算精度较高且具有自启动的特点。对于用龙格库塔法求解单个常微分方程和求解常微分方程组的思路基本相
2、似(注意一点一个微分方程组是常微分方程组即表明微分方程中的各阶导数都是对同一个变量求导,例如可以把各个量对时间求导得到一个常微分方程组,如果一个微分方程组中的有对不同变量的导数那么这个方程组就成了偏微分方程组),都是根据泰勒展开得到其迭代计算形式,其基本思想都是按照一定原则求取当前点附近一些点的斜率,通过这些斜率的线性组合作为当前点处的斜率,进行递推求解。2、数学模型2.1常微分方程初值问题的数学模型(1)式中为的已知函数,为给定的初始值。常微分方程的数值解法的任务就是要求出函数值在微分自变量取如下序列:时的值的近似值,一般情况下都采取等距结点的方式,即:,其中为相邻两结点的距离,称为步长。2
3、.2常微分方程组初值问题的数学模型(2)式中为的已知函数,为初始值,利用数值解法求解常微分方程组的任务就是求出函数值在微分自变量取如下序列:时其值的近似值3龙格库塔法的迭代形式龙格库塔法的迭代形式推导在相关的数值分析书籍中均有详细的讲解,这里就不进行推导直接给出常用的四阶龙格库塔法的形式作为例子。3.1单个常微分方程的迭代形式对应于2.1节的数学模型有:(3)通过上面的迭代形式可知实际上就是如下四个点(点1)(点2)(点3)(点4)处函数关于的斜率。这四个点的确定又是一环扣一环的,比如要求得点2的位置,就必须先求得在点1处的斜率3.2常微分方程组的迭代形式对应于2.2节的数学模型有:(4)对比
4、于单个微分方程的龙格库塔迭代形式,可以看出整个求解的过程完全一样,只不过在每步算斜率的时候需要一次性的解算完(方程的个数)个,即是要把的斜率一次性算完,才能继续往下递推运算。在求解常微分方程组初值问题的情况下,可将常微分方程组写成向量的形式,这样更加便于理解,下面用图形来展示整个求解过程:图1龙格库塔法的求解过程沿着箭头的方向依次解出各个向量,即可依次递推下去,下面的表示为点处的斜率,处情形类似。4算例4.1问题描述利用四阶龙格库塔法求解初值问题取步长为。4.2问题分析为了和第3节的数学模型对应,我们令,,则上面的问题就变为和前面的模型一样的形式:4.3C程序代码#include5、h>#include#defineM2//维数,方程组中方程的个数//算斜率的右端函数voidRightSlopeFun(double*pRightK,//算出来的右端斜率项[out]doubleiX,//计算点的微分自变量值double*iY//计算点的函数初始值){pRightK[0]=iX*iY[0]-iY[1];pRightK[1]=(iX+iY[0])/iY[1];}//龙格库塔法解算程序doubleLgktSoulution(double*pOY,//返回函数值[out]doubleiX,//积分自变量输入值double*iY,//初始值doubleh,//单步步长6、intCount=1//解算步数,默认形参为1,调用时不给此//参数赋值,则只解算一步){doubleK1[M];//第一个斜率doubleK2[M];//第二个斜率doubleK3[M];//第三个斜率doubleK4[M];//第四个斜率doubletempY[M];//保存Y的中间值的临时空间inti;//计数器doubletempX;//临时保存x的空间tempX=iX;//载入初始的x值for(i=0;i0){RightSlopeFun(K1,tempX,tempY)7、;//解算斜率K1tempX+=h/2;for(i=0;i
5、h>#include#defineM2//维数,方程组中方程的个数//算斜率的右端函数voidRightSlopeFun(double*pRightK,//算出来的右端斜率项[out]doubleiX,//计算点的微分自变量值double*iY//计算点的函数初始值){pRightK[0]=iX*iY[0]-iY[1];pRightK[1]=(iX+iY[0])/iY[1];}//龙格库塔法解算程序doubleLgktSoulution(double*pOY,//返回函数值[out]doubleiX,//积分自变量输入值double*iY,//初始值doubleh,//单步步长
6、intCount=1//解算步数,默认形参为1,调用时不给此//参数赋值,则只解算一步){doubleK1[M];//第一个斜率doubleK2[M];//第二个斜率doubleK3[M];//第三个斜率doubleK4[M];//第四个斜率doubletempY[M];//保存Y的中间值的临时空间inti;//计数器doubletempX;//临时保存x的空间tempX=iX;//载入初始的x值for(i=0;i0){RightSlopeFun(K1,tempX,tempY)
7、;//解算斜率K1tempX+=h/2;for(i=0;i
此文档下载收益归作者所有