一阶常微分方程数值解的C语言编程实现.doc

一阶常微分方程数值解的C语言编程实现现实问题与分析学科自身的发展，使得微积分及其涉及领域内出现了寻求数值解的诸多问题(在其他数学领域亦是如此)。这是由于多种原因造成的，比如现实问题中函数表达式往往并不存在，即使可以拟合也并没有得到表达式的必要且表达式的讨论可能出现较大误差。再比如，很多复杂函数或函数约束条件很难使用甚至不能使用一般的分析工具，像微分方程存在不可积、不可解的大量实例。所以数值解的相关理论不断产生、完善。特别是现代以来，计算机的迅速发展为解决数值解问题提供了算法的计算机程序实现可能，并提供了更多的求解思路。Mathematica，MATLAB，Maple是最常用的数学软件，对应于Mathematica语言，MATLAB语言，Maple语言。当然，一般使用C语言也可以比较便捷地解决很多分析学科问题的数值解。举一元函数微积分学一些熟知的数值解例子：数值导数常见有图解微分法、差商公式、插值多项式、拉格朗日公式(拉格朗日插值公式得来)、马尔科夫公式(牛顿插值公式得来)、等距公式、三次样条函数等等，数值积分有内插求积公式(最为基础的有梯形法、辛普森公式(抛物线法)、高斯型求积公式等等。这里，考察一阶导数已解出的一阶常微分方程的初值问题的数值解。问题给出：常微分方程数值解问题中，一阶常微分方程与一阶常微分方程组(一阶导数均解出形式)的数值解法是最基础的，也是非常常用的。一般数值解法有欧拉方法及改进的欧拉方法、龙格-库塔方法、阿达姆斯方法。现给出前两种方法，它们都是单步法，而阿达姆斯法是线性多步法，各有优缺点(此外计算机上常会加以使用浮动步长法，如自动变步长计算方法)。单步法的通用表达式是拉格朗日中值定理的形式(差分代替微分是一种重要的近似思想，反过来用微分代替差分亦重要)：欧拉方法：改进的欧拉方法(预报校正法)：龙格-库塔方法：二阶龙格-库塔方法：四阶龙格-库塔方法：容易看出，欧拉方法及改进的欧拉方法都是特殊的龙格-库塔方法。对于龙格-库塔公式，常数的确立是通过单、双变量泰勒级数展式推理得到的，不再叙述。另外，微分方程数值解的相容性条件、收敛性问题、精度判断、误差分析等内容也不叙述。只说明龙格-库塔方法通常使用上述的二阶、四阶形式。下面给出一例：这是变量可分离方程，容易积分得到：现写出该一阶常微分方程数值解法的一个C语言程序，输出结果是分别使用欧拉方法、预报校正法、二阶龙格-库塔法、四阶龙格-库塔法得到的函数近似值组成的表格打印。程序如下：#include stdio.h double f(double x,double y)return y*(1-y*y)；void main()int n,i；float x,ye,yeu,yeul,yr,k11,k21,yru,k1,k2,k3,k4,x0,y0,h；printf(n请输入初值x=)；scanf(%f,&x0)；printf(y=)；scanf(%f,&y0)；printf(步长h=)；scanf(%f,&h)；printf(总步数n=)；scanf(%d,&n)；x=x0；ye=yeu=yr=yru=y0；printf(n%7s欧拉方法改进的欧拉方法二阶龙格-库塔方法四阶龙格-库塔方法n,x)；for(i=1；i=n；i+)ye+=h*f(x,ye)；yeul=yeu+h*f(x,yeu)；yeu+=(h/2)*(f(x,yeu)+f(x,yeul)；k11=f(x,yr)；k21=f(x+h/2,yr+(h/2)*k11)；yr+=h*k21；k1=f(x,yru)；k2=f(x+h/2,yru+(h/2)*k1)；k3=f(x+h/2,yru+(h/2)*k2)；k4=f(x+h,yru+h*k3)；yru+=(h/6)*(k1+2*k2+2*k3+k4)；x+=h；printf(%8.1f%13.4f%13.4f%17.4f%19.4fn,x,ye,yeu,yr,yru)；printf(nn)；执行结果如下：这个程序涉及语句非常简单，而且具有很强的可移植性，即改变第四行return语句规定的被调函数返回值，代以其它一阶常微分方程已解出的导数表达式(注意相关变量的数据类型声明以及特殊常量的定义，此外使用数学函数还需调用头文件math.h)，此程序依旧可以执行上述四种数值解法，解决代换的方程的函数值求近似解问题。要说明的是，这里步长应精确到十分位选取，如果要进一步精确，可在倒数第四行%8.1f处将小数位数指定由1变成所需位数，但不可过大而扰乱表格排列。本例属于基本类型的一阶微分方程，直接积分已经得到结果。将精确解编入程序，编写打印数值解误差的程序。突出改变重点地列出新程序语句：#include math.h#include stdio.h double f(double x,double y)return y*(1-y*y)；double Y(double x)return sqrt(4*exp(2*x)/(4*exp(2*x)-3)；void main()float y；x+=h；y=Y(x)；printf(%8.2f%13.4f%13.4f%17.4f%19.4fn,x,fabs(ye-y),fabs(yeu-y),fabs(yr-y),fabs(yru-y)；执行如下：可以看出四阶龙格-库塔公式整体误差始终最小，它有四阶精度。而有二阶精度的改进的欧拉方法、二阶龙格-库塔法误差相似，一阶精度的欧拉方法误差较大，但也可大致估计函数性态。为了综合函数近似值、整体误差分析，对程序作以语法修改，得出最后给出的程序。对于一般的一阶导数已解出一阶常微分方程，此程序均可以使用，且调整一阶导数表达式、原函数表达式(研究数值解误差分析时需要输入原函数表达式；如果只是为了求得方程数值解，表达式处随意输入符合定义的返回值，更改for(s=1；s=2；s+)循环控制条件为s 2)，只需变动前几行中相应的两个return语句；而步长精确度修改方法同上。程序如下：#include math.h#include stdio.h double f(double x,double y)return y*(1-y*y)；double Y(double x)return sqrt(4*exp(2*x)/(4*exp(2*x)-3)；void main()int n,s,i；float x,y,ye,yeu,yeul,yr,k11,k21,yru,k1,k2,k3,k4,x0,y0,h；printf(n请输入初值x=)；scanf(%f,&x0)；printf(y=)；scanf(%f,&y0)；printf(步长h=)；scanf(%f,&h)；printf(总步数n=)；scanf(%d,&n)；for(s=1；s=2；s+)if(s 2)printf(n函数近似值：n)；else printf(整体误差：n)；printf(n%7s欧拉方法改进的欧拉方法二阶龙格-库塔方法四阶龙格-库塔方法n,x)；x=x0；ye=yeu=yr=yru=y0；for(i=1；i=n；i+)ye+=h*f(x,ye)；yeul=yeu+h*f(x,yeu)；yeu+=(h/2)*(f(x,yeu)+f(x,yeul)；k11=f(x,yr)；k21=f(x+h/2,yr+(h/2)*k11)；yr+=h*k21；k1=f(x,yru)；k2=f(x+h/2,yru+(h/2)*k1)；k3=f(x+h/2,yru+(h/2)*k2)；k4=f(x+h,yru+h*k3)；yru+=(h/6)*(k1+2*k2+2*k3+k4)；x+=h；y=Y(x)；if(s 2)printf(%8.1f%13.4f%13.4f%17.4f%19.4fn,x,ye,yeu,yr,yru)；else printf(%8.1f%13.4f%13.4f%17.4f%19.4fn,x,fabs(ye-y),fabs(yeu-y),fabs(yr-y),fabs(yru-y)；printf(nn)；执行情况：现通过一张表格，说明满足初值条件y(0)=2的y=y(x)在x+时y1，且y在x是很小的正数(最高位为个位)时就已经可以开始高度精准地当作1-表格中初值为y(5)=1即可代替y(0)=2。实质上，从既得解析式可以迅速观察得到，x+时函数以1为极限，且函数是从y极大快速降至持续贴近1，而这种下降是在x为很小的正数(最高位为个位)时。为了形象地表达其趋于1的情形，附带给出该方程积分曲线。待求解函数的极限性质便清晰得出。最后，再随