数值计算一阶常微分方程求解实现

原创作品，出自 “晓风残月xj” 博客，欢迎转载，转载时请务必注明出处（http://blog.csdn.net/xiaofengcanyuexj）。

由于各种原因，可能存在诸多不足，欢迎斧正！

二、引言

常微分方程是在解决实际问题的过程中产生的，微分方程的研究又促进实际问题的解决，从这可以看出微分方程的重要性。关于微分方程的问题一般比较复杂，我们在这里只选取其中较为简单的一阶常微分方程进行研究与讨论。

很多情况下实际问题是能通过抽象建模转化为一阶常微分方程形式的，然而即便给出了一阶常微分方程，想通过它来求解函数在某个点的值也并不是一件容易的事。因为即使给出了一阶微分方程及初值条件，即便不考虑研究者的知识水平，往往也很难用现有数学理论和方法来导出原函数，而数值计算则为解决一阶常微分方程问题提供了强有力的理论依据。我们知道，计算机处理的问题具有离散、不连续的基本特征，与一阶常微分方程的数值计算问题高度吻合，可以借助计算机自动、高速运算的特性来有效解决一阶常微分方程数值分析过程中繁琐的计算问题。

因此，有必要研究运用计算机解决带初始条件的一阶常微分方程数值计算问题。

三、算法描述

1，Euler method

公式：Yn+1=Yn+h*f(tn,Yn)是一个递推公式，由初始条件Y0=C,即可计算Y1的近似值

Y1=Y0+ h*f(t0,Y0)

Y2=Y1+h*f(t1,Y1)

……………………

Yn+1=Yn+h*f(tn,Yn)

不断的进行迭代，即可得到Yn+1的值。

显然可以得递推的计算公式

Y0=C

Yn+1=Yn+h*f(tn,Yn) n=0,1,2,3……

图示曲线是给定的y=y(x),欧拉方法就是用图中的折线近似地代替曲线y=y(x)，因此欧拉法有时也称折线法。

Euler method我们的具体实现：

inc=(t-start_t)*1.0/(n*1.0);

若choice等于1(注释：start_t=0,t:用户输入的数字，n=pow(2.0,j*1.0),j初始值为1，最大值14)，则使用已经初始化微分方程

dy/dt=t*t*y,y[0]=1进行计算

由公式y[tn+1]=y[tn]+inc*f（tn,y[tn]）

（1）通过迭代求得y(n),存入sum[j],若在迭代过程中求的f(ti,y[ti])趋近无穷时则表明在 ti 不可导,退出循环

（2） j=j+1，执行步骤（1），如果sum[j]-sum[j-1]< infinitesimal(注释：一个极小的数)或者j>15,则退出循环，sum[j]即为微分方程在t点的近似解

若choice等于2，则选择用户自己输入的微分方程进行计算数值，为了简便起见，n= t/INC+1 (INC是一个比较小的常数) ，j设置为14，即执行一次循环，即可得出在t点的数值。

2,Blackward Euler

导数Yn’近似值表达式为：

Yn’≈（Yn-Yn-1）/h

由此得到方程

Yn+1=Yn+h*f(tn+1,Yn+1)

图形化的显示如下：

Y=y(x)为途中的曲线，折线是由向后欧拉法得到的

Blackward Euler我们的具体实现：

inc=(t-start_t)*1.0/(n*1.0);

若choice等于1(注释：start_t=0,t:用户输入的数字，n=pow(2.0,j*1.0),j初始值为1，最大值14)，则使用已经初始化微分方程

dy/dt=t*t*y,y[0]=1

进行计算

由公式y[tn+1]=y[tn]+inc*f（tn+1,y[tn+1]）

令f[y([tn], tn, y[tn+1])=y[tn+1]-y[tn]-inc*f（tn+1,y[tn+1]）

（1）调用二插法函数Bisection1()求方程y[n+1]=y[n]+inc*f（t[n+1],y[n+1]）的解y[n+1],若y[i+1]趋近无穷时则表明在 ti点不可导,退出循环，否则存入sum[j]

（2）j=j+1，执行步骤（1），如果sum[j]-sum[j-1]< infinitesimal(注释：一个极小的数)或者j>15,则退出循环，sum[j]即为微分方程在t点的近似解

3,Crank-Nicholson

如果我们使用中间差异而不是向前差的欧拉法或后向差分的向后欧拉法,我们会得到一个二阶方法

Yn+1/2’≈(Yn+1-Yn)/h

替换成微分方程为：

（Yn+1-Yn）/h≈f(tn+1/2,Yn+1/2)

由此得到

（Yn+1-Yn）=1/2*[f(tn+1,Yn+1)+f(tn,Yn)]*h

Crank-Nicholson我们的具体实现：

inc=(t-start_t)*1.0/(n*1.0);

若choice等于1(注释：start_t=0,t:用户输入的数字，n=pow(2.0,j*1.0),j初始值为1，最大值14)，则使用已经初始化微分方程

dy/dt=t*t*y,y(0)=1

进行计算

由公式y[tn+1]=y[tn]+0.5*inc*（f（tn,yn）+f（tn+1,y[tn+1]））

令f（y[tn], tn, y[tn+1]）=y[tn+1]-y[tn]- 0.5*inc*[f（tn,yn）+

f(tn+1,y[tn+1])]

（1）调用二插法函数Bisection1()求方程y[n+1]=y[n]+1/2*inc*[f(t[n+1],y[n+1])+f(t[n],y[n])]的解y[n+1],,若y[i+1]趋近无穷则表明在 ti点不可导,退出循环, 否则存入sum[j]

（2） j=j+1，执行步骤（1），如果sum[j]-sum[j-1]< infinitesimal(注释：一个极小的数)或者j>15,则退出循环，sum[j]即为微分方程在t点的近似解

4,Multistep methods

作为一种替代方法,其精度逼近的导数可以通过使用一个线性组合的附加点改善。如果s=2有

Yn’=fn

Yn’’=∽(fn-fn-1)/h

所以我们可以构造一个二阶泰勒级数展开的方法是

Yn+1=Yn+h*Yn’+1/2*h*h*Yn’’=Yn+1/2*h*(3*fn-fn-1)

对于s = 3我们也有Y”n,所以可以利用一个二阶片面的有限差分来近似Y“n = f’=(3 fn -4 fn-1 + f n -2)/ 2h和包括三阶Y”“n = f“n =(fn - 2 fn-1 + fn-2)/(h*h)获得

Yn+1=Yn+h*Yn’+1/2*h*h*Yn’’+1/6*h*h*h*Yn’’

Yn+1=Yn+1.0/12*(23*fn-16*fn-1+5*fn-2)

Multistep methods我们的具体实现：

Inc1=(t-start_t)*1.0/(n*1.0); 定义derivative[]数组存放导数值

若choice等于1(注释：start_t=0,t:用户输入的数字，n=pow(2.0,j*1.0),j初始值为1，最大值14)，则使用已经初始化微分方程

dy/dt=t*t*y,y(0)=1

进行计算

由公式y(tn+1)=y(tn)+(1.0/12)*inc*(23*fn-16*fn-1+5*fn-2)

知道需要先求处derivative[0], derivative[1], derivative[2]的值，显然derivative[0] =0,由Euler method可以得到root[0]=Euler_method(start_t),temp= start_t +1*inc1,derivative[1] =f(temp,root[0]);temp=start_t+2*inc1,同理可以得出derivative[2]的值。

以下用迭代法求数值定义变量i=3

（1）j=1

（2）temp=start_t+inc1*i;

after+=1.0/12.0*inc1*(23*derivative[2]-16*derivative[1]+5*derivative[0]);

derivative[0]=derivative[1];

derivative[1]=derivative[2];

tempderivative=example_equation(temp,le);

若导数tempderivative过大或过小，就认为该点不可导，返回

（3） i=i+1,继续执行（1），如果i>n,则退出循环，after存入

sum[j]

（4） j=j+1,若j>=15或者sum[j]-sum[j-1] <infinitesimal（注释：infinitesimal：一个极小的数），则退出循环，sum[j]即为微分方程在t点的近似解，否则转到（2）继续执行

若choice等于2，则选择用户自己输入的微分方程进行计算数值，为了简便起见，n= t/INC+1 (INC是一个比较小的常数表示步====) ，j设置为14，即执行一次循环，即可得出在t点的数值。

以下图示为用三种方法得到的折线，以及原函数曲线。

四、算法分析与运行结果分析

1.Euler method

Euler method的误差为O(Dt2)，算法思想简单，易于理解，同时实现起来也比较容易，以下是我们组的运行结果（以曾老师所给实例为准）

2.Backward Euler的误差为O(Dt2)，算法思想也简单，也易于理解，同时实现起来也比较容易，以下是我们组的运行结果（以曾老师所给实例为准）

3.Crank-Nicholson的误差O(Dt2)，不同之处在于导数的表示方法不一样，Crank-Nicholson是用小区间中点的导数来表示的，而且有个技巧近似的用两端点导数的算术平均值来估计中点导数，有效减少运算量

4.Multistep methods以s=3展开，误差为O(Dt3)

五、我们的特色

1.菜单的设置

本程序设计了菜单界面，可以选择使用1、2、3号功能。

1号功能代表直接对∫dy/dt=t*t*y其中y(0)=1的一阶常微分方程进行数值计算，我们有意识的把初始值t0到所求自变量t的区间由小到大的分为不同14种，这样可以验证分的越小，越接近真实值这一基本原则，而却由于函数简单且可分离，我们求出了原函数，通过原函数求出精确值，然后进行对比，验证算法的正确性，分析误差。

2号功能则是对用户输入的函数进行数算。此时不知道精确值，无法像1号功能那样进行验证性计算，只能求出近似值，在有1功能的前提下是可以保证算法实现程序的正确性与可靠性的。此外，1号功能直接计算dy/dt=t*t*y其中y(0)=1的一阶常微分方程，资源占用小，运算速度快；而2号功能用户输入的函数则要通过繁琐的字符串识别与表达式的计算,涉及进出栈需占用大量系统资源，时间和空间尤其是时间资源消耗较大，所以我们最后直接计算出结果，不把初始值t0到所求自变量t的区间由小到大的划分过程展示出来。

此外，可以通过2功能输入dy/dt=t*t*y其中y(0)=1的一阶常微分方程与1功能相互验证，说明程序的正确性与健壮性。如下：

2.实现了用户输入函数

我们的程序可以根据用户的具体计算需要输入函数（由于写程序初期有些问题未意识到，所以并不能覆盖所有的函数，但较之在代码中封装单一dy/dt=t*t*y其中y(0)=1函数有很大改善），然后通过自己写的头文件"识别函数程序4.h"，完成函数的识别工作。具体工作为以字符串形式输入用户函数，然后结合表达式的转换与计算知识将输入的中缀表达式转换为后缀表达式存储在自定义堆栈中，在此基础上根据用户输入的积分区间,调用<math>库中的函数计算后缀表达式的值（具体实现见源代码）。

我们组可以识别的函数举例（有些函数可能不常见，这里不做解释）

fabs cos sin tan asin acos atan

sinh cosh tanh exp floor log

以及它们与x的任意组和，复杂多项式等等，在注意了输入格式的前提下基本可以实现绝大多数函数的带初始条件的一阶常微分方程数值计算。

下面是一组用户随机函数运用课本带初始条件的一阶常微分方程数值计算的4种思想得到的结果，有兴趣的同学可以验证一下：

用户函数1:

dy/dt=(cos(x)-y)/t其中y(PI)=1

运算结果如下：

3.将计算结果与实际对照

如用户函数2:

dy/dt=(2*t*y-y)/t^2其中y(1)=e

下面是经过手动数值计算得到的原函数表达式，然后代自变量值通过计算器精确计算所得

下面是我们的程序计算所得，

可以看出存在一定误差，但都是可以接受的。

通过这种数值计算与手动计算结合的方式可以增加说服力，证明程序算法实现的正确性。

4.Backward Euler与Crank-Nicholson的处理

Backward Euler 要处理

Crank-Nicholson 要处理

这两种方法都要解决方程求根问题，即在其他条件已知的情况下求y[n+1],我们组是学以致用采用前面关于方程求根中的二分法Bisection，之所以选二分法Bisection，是因为二分法Bisection简单可靠，易于实现而容易理解，可以将上面两个表达式依次转化为：

y[n+1] -（y[n]+inc*f(t[n+1],y[n+1])）=0

y[n+1] -（y[n]+1/2*inc*[f(t[n+1],y[n+1])+f(t[n],y[n])）=0

然后通过二分法Bisection求y[n+1]。

此种方法虽然有它的先天不足（下面介绍），但可人为控制以保证程序的正确性与健壮性。

5.Multistep methods的处理

此方法要知道f[n+1] ，f[n] ，f[n-1],所输入方程的右边就导数，求出倒数这不是问题，但计算只能从第3个小区间点开始依次向后推进，前面的f[0] ，f[1] ，f[2]未知，需要程序求出。，可以先通过Euler method求得y[1],y[2](其中y[0]为初始条件不用求)，然后代人一阶常微分方程的右边（此处由我们的输入形式决定，默认左边只是dy/dt的形式，右边即为导数）求得导数f[0] ，f[1] ，f[2]。

但这样的方法算出的结果有误差，精确度与效率不如Crank-Nicholson。

6.误差分析

7.不可导函数的处理

一阶常微分方程数值不同于前面的积分的数值计算问题，前者相对复杂，计算中存在大量不确定因素，一阶常微分方程的原函数可能不可导（不连续，左右导数不相同等等），如不考虑此类问题在数值计算时是可能引发灾难性错误的（即在基本计算方法正确的情况下，不知不觉的计算结果完全错误）

如输入用户函数3：

dy/dt=y^2其中y(1)=2

其原函数图像如下：

可以看出原函数在替t=1.5不可导；这样如果不做任何处理必然导致错误

经处理我们组的计算结果如下：

可以看出：

Euler method提示在t=1.50987处不可导；

Backward Euler method，Crank Nicholsom显示结果为1e100,，这是我们组设置的二分法Bisection右边区间初始值。我们将二分法Bisection左右区间初始值分别设置为NEGATIVE=-1e+50，NO_NEGATIVE=1e+50，因为二分法Bisection求方程根时要给定初始区间范围；但一阶常微分方程数值计算区间是动态的，随着t的变化而变化，而且可能突然就变得非常大，如原函数为指数形式，因此依据前一步的计算结果设置二分法Bisection的初始区间是不可取的，所以有必要把二分法Bisection的初始区间设置得非常大，但这样会导致计算机时间和空间资源的消耗，但为了保证正确性与健壮性有必要将二分法Bisection的初始区间设置得非常大，在正确性与效率产生矛盾时我们选择了牺牲效率保证正确性；

Multistep methods提示在t=1.5196201处不可导；

这些提示信息的输出有效的提醒了用户当前程序计算结果所处的状态。

说明：

以上所选一阶常微分方程都是可以找到原函数的，这样做是为了将计算机数值计算结果与手动计算结果比较，验证算法实现程序的正确性与健壮性。也可以输入不易找到原函数的一阶常微分方程，但计算结果的正确无法得到确认。

六、总结

这次程序我们有一个问题由于各种原因未能完全解决。如下：

Backward Euler

Crank-Nicholson

对于以上两个式子的处理，由于实现了用户输入函数，所以我们没有直接分离求出Y[n+1]的表达式（这点是比较繁琐的），而是选择二分法Bisection，撇下效率不说，正确性上也有问题。如果原函数有多个根，通常二分法Bisection在正负无穷区间上只能求得其中一个根，这个根不一定就是所要求的Y[n+1]。关于这点我们会继续改进。

通过这学期计算方法的学习与交流，我们加深了对计算机解决实际问题方法的认识：

1.提出问题

2.建立模型

3.理论分析（这点很重要，就拿这次而言，这4种带初始条件的一阶常微分方程数值计算的基本思想，为带初始条件的一阶常微分方程数值计算的实现提供了理论依据，同时也证明了其方法的正确性与可行性，使我们编程时有了明确的指导思想不至于盲干）

4.编程实现

#include<iostream>//头文件
#include<math.h>
#include<string>
#include<stack>
using namespace std;char expression[200];bool operator_priority(char temp_operator,char expression)
//若栈顶运算符temp_operator优先级比当前表达式的高，则返回true,否则返回false
{if(temp_operator=='-'||temp_operator=='+'){if(expression=='-'||expression=='+'||expression==')'||expression=='#')return true;else if(expression=='/'||expression=='*'||expression=='('||expression=='^')return false;}else if(temp_operator=='/'||temp_operator=='*'){if(expression=='-'||expression=='+'||expression=='/'||expression=='*'||expression==')'||expression=='#')return true;else if(expression=='('||expression=='^')return false;}else if(temp_operator=='('){if(expression!=')')return false;else true;}else if(temp_operator==')')return true;else if(temp_operator=='#')return false;else if(temp_operator=='^')/{if(expression=='(')return false;else true;}
}bool is_operator(char expression)
//判断expression是否为运算符，若是返回true，否则返回false
{if(expression=='-'||expression=='+'||expression=='*'||expression=='/'||expression==')'||expression=='('||expression=='^'||expression=='#')return true;else return false;
}int len;
string post;
string str[200];
void mid_trans_post()
//将中缀表达式转化为后缀表达式存储在str[200]
{int i;gets(expression);stack<char>Stack_operator;Stack_operator.push('#');i=0;//  cout<<"expression="<<expression<<endl;while(i<strlen(expression)){if(expression[i]>='0'&&expression[i]<='9')//代表的是数据post=post+expression[i];else if(expression[i]>='a'&&expression[i]<='z')//代表的是数据//post=post+expression[i];else if(i>0&&!is_operator(expression[i-1])&&is_operator(expression[i])){str[len++]=post;/*cout<<" post="<<post;*/post.erase(0,post.size());}if(is_operator(expression[i]))//代表是操作符{
leap:   char temp_operator=Stack_operator.top();        if(temp_operator=='('&&expression[i]==')'){   Stack_operator.pop();}else if(temp_operator=='#'&&expression[i]=='#'){Stack_operator.pop();break;}else if(temp_operator==')'&&expression[i]=='('){   cout<<"你的表达是有错!"<<endl;exit(0);}else if(!(operator_priority(temp_operator,expression[i]))){Stack_operator.push(expression[i]);}else if(operator_priority(temp_operator,expression[i])){   Stack_operator.pop();str[len++]=temp_operator;while(!Stack_operator.empty()){goto leap;}}}i++;}
}double post_Exp(double iteration_x,double iteration_y)
//计算代值iteration_x之后后缀表达式的值
{double x1,x,x2;int i;stack<double>Stack_operand;for(i=0;i<len;i++){x=0;if(str[i][0]>='0'&&str[i][0]<='9'){int j=0;while(j<str[i].size()){x=x*10+str[i][j]-'0';j++;}// cout<<" x="<<x<<" ";Stack_operand.push(x);}else if(str[i][0]>='a'&&str[i][0]<='z'){double temp;if(str[i]=="x")temp=iteration_x;else if(str[i]=="y")temp=iteration_y;else if(str[i]=="fabs")temp=fabs(iteration_x);else if(str[i]=="cos")temp=cos(iteration_x);else if(str[i]=="sin")temp=sin(iteration_x);else if(str[i]=="tan")temp=tan(iteration_x);else if(str[i]=="asin")temp=asin(iteration_x);else if(str[i]=="acos")temp=acos(iteration_x);else if(str[i]=="atan")temp=atan(iteration_x);//  else if(str[i]=="atan2")temp=atan2(1,2);else if(str[i]=="sinh")temp=sinh(iteration_x);else if(str[i]=="cosh")temp=cosh(iteration_x);else if(str[i]=="tanh")temp=tanh(iteration_x);else if(str[i]=="exp")temp=exp(iteration_x);else if(str[i]=="floor")temp=floor(iteration_x);//  else if(str[i]=="fmod")temp=fmod(5,3);else if(str[i]=="log")temp=log(iteration_x);//  else if(str[i]=="log10")temp=log10(105);Stack_operand.push(temp);}else {x2=Stack_operand.top();//   cout<<" x2="<<x2<<" ";Stack_operand.pop();x1=Stack_operand.top();//  cout<<" x1="<<x1<<" ";Stack_operand.pop();switch(str[i][0]){case '+':x1+=x2;break;case '-':x1-=x2;break;case '*':x1*=x2;break;case '/':x1/=x2;break;case '^':x1=pow(x1,x2);}Stack_operand.push(x1);}}return Stack_operand.top();
}void menu()
//输出提示并调用相关函数进行计算
{int i;cout<<"*****************************************************************"<<endl;cout<<"Please input the  function(care for the format):"<<endl;len=0;mid_trans_post();cout<<"followings is the infix of the  function:"<<endl;for(i=0;i<strlen(expression)-1;i++){if(expression[i]!=' ')cout<<expression[i];}cout<<endl<<"followings is the suffix of the  function:"<<endl;for(i=0;i<len;i++){cout<<str[i]<<" ";}cout<<endl<<"*****************************************************************"<<endl;}
/*
int main()
{while(1){menu();getchar();}return 0;
}*/

#include<iostream>
#include"识别函数程序4.h"
#include<iomanip>
#include<math.h>
#define MAX 100
#define  infinitesimal 1e-5//表示两数相差很小
#define NEGATIVE -1e-50///表示负无穷大
#define NO_NEGATIVE -1e+50///表示正无穷大
using namespace std;int sum_num=0;/记录数组sum[]的长度
int choice;
double sum[MAX];
double inc;///表示小区间长度double start_t=0;
double start_y=1;///表示初始条件t0,y0//测试例子
//          dy/dt=t^2*y
//          y(0)=1void print(int num)
//以固定格式输出
{int i;cout<<"          index    number of interval                root  "<<endl;for(i=1;i<num;i++){cout<<setw(15)<<i<<setw(22)<<(int)pow(2.0,i)<<setw(20)<<setprecision(8)<<sum[i]<<endl;      }
}double example_equation(double t,double y)
//测试例子的倒数函数（右边）
{return t*t*y;
}Euler method法求一阶微分方程//
double Euler_method(double t)
//Euler method法求一阶微分方程
{int j;double temp;double before,after;double n;before=start_y;反复进行直到前后两次绝对值之差小于infinitesimal或超过15次退出for(j=1;j<15;j++){n=pow(2.0,j*1.0);after=start_y; inc=(t-start_t)*1.0/(n*1.0);double i=1;//通过迭代求得Y[n]while(i<=n){temp=start_t+inc*i;      if(1==choice)after=after+inc*example_equation(temp-inc,after);else if(2==choice)after=after+inc*post_Exp(temp-inc,after);       i++;}     sum[j]=after;///找到近似解if(j>1&&fabs(after-before)<infinitesimal){j++;break;}before=after;}sum_num=j;return after;
}Backward Euler method法求一阶微分方程//
double f(double before,double t,double y)
//求函数y[n+1]-(y[n]+inc*f(t[n+1],y[n+1]))的值
{if(1==choice)return y-(before+inc*example_equation(t,y));else return y-(before+inc*post_Exp(t,y));
}double Bisection(double before,double t,double left,double right)
//二插法求方程y[n+1]=y[n]+inc*f(t[n+1],y[n+1])的解y[n+1],
{double mid;while(fabs(right-left)>infinitesimal){mid=(left+right)/2;if(f(before,t,left)*f(before,t,mid)>=0) left=mid;else if(f(before,t,right)*f(before,t,mid)>=0)right=mid;        }return mid;
}double Backward_Euler_method(double t)
//Backward Euler method法求一阶微分方程
{int j;double temp;double before,after,n;before=start_y;反复进行直到前后两次绝对值之差小于infinitesimal或超过15次退出for(j=1;j<15;j++){n=pow(2.0,j*1.0);after=start_y;    inc=(t-start_t)*1.0/(n*1.0);double i=1;//通过迭代求得Y[n]while(i<=n){temp=start_t+inc*i;after=Bisection(after,temp,NEGATIVE,NO_NEGATIVE);i++;} sum[j]=after;///找到近似解if(j>1&&fabs(after-before)<infinitesimal){j++;break;}before=after;}  sum_num=j; return after;
}Multistep methods法求一阶微分方程///
double derivative[3];
double Multistep_methods(double t)
//Multistep methods程法求一阶微分方程
{int j;double temp;/double before,after,n;double inc1;//    int step;//cout<<"请输入步数:"///cin>>step;before=start_y;反复进行直到前后两次绝对值之差小于infinitesimal或超过15次退出for(j=1;j<15;j++){n=pow(2.0,j*1.0);after=start_y;inc1=(t-start_t)*1.0/(n*1.0);///依次求得f[0],f[1],f[2]double root[3];temp=start_t; derivative[0]=0;root[0]=Euler_method(temp);temp=start_t+inc1;root[1]=Euler_method(temp);if(1==choice)derivative[1]=example_equation(temp,root[0]);else if(2==choice)derivative[1]=post_Exp(temp,root[0]);temp=start_t+2*inc1;root[2]=Euler_method(temp);if(1==choice)derivative[2]=example_equation(temp,root[1]);else if(2==choice)derivative[2]=post_Exp(temp,root[1]);after=root[2];double i=3;//通过迭代求得Y[n]while(i<=n){temp=start_t+inc1*i;double le=after;//利用公式y[n+1]=y[n]+1/12*inc*(23*f[n]-16*f[n-1]+5*f[n-2])after+=1.0/12.0*inc1*(23*derivative[2]-16*derivative[1]+5*derivative[0]);/求运算后的f[n],f[n-1],f[n-2]derivative[0]=derivative[1];derivative[1]=derivative[2];if(1==choice)derivative[2]=example_equation(temp,le);else if(2==choice)derivative[2]=post_Exp(temp,le);i++;}    sum[j]=after;///找到近似解if(j>1&&fabs(after-before)<infinitesimal){j++;break;}before=after;}  sum_num=j;return after;
}//主函数
int main()
{double t;while(1){int tag=1;///根据菜单提示选择相应功能///int flag=0;do{if(flag)cout<<"Your choice is illegal,please choose again!"<<endl;cout<<"==============================menu=================================="<<endl;cout<<"              1.Integral verification on the books"<<endl;cout<<"              2.Input new function and do numerical integration"<<endl;cout<<"              3.Quits the program"<<endl;cout<<"Please choose [ ]\b\b";cin>>choice;getchar();cout<<"===================================================================="<<endl<<endl;flag++;}while(choice!=1&&choice!=2&&choice!=3);    if(2==choice) menu();else if(3==choice)exit(0);cout<<"Please input variable t:";cin>>t;                   cout<<"followings is the root of the  function:"<<endl<<endl;double root;cout<<"******************************************************************"<<endl;cout<<"              "<<tag++<<".Here is the result of Euler method"<<endl;root=Euler_method(t);print(sum_num);cout<<"Here is the root of Euler method:"<<setw(20)<<setprecision(8)<<root<<endl<<endl;cout<<"******************************************************************"<<endl;cout<<"              "<<tag++<<".Here is the result of Backward Euler method"<<endl;root=Backward_Euler_method(t);print(sum_num);cout<<"Here is the root of Backward Euler method:"<<setw(20)<<setprecision(8)<<root<<endl<<endl;cout<<"******************************************************************"<<endl;cout<<"              "<<tag++<<".Here is the result of Multistep methods"<<endl<<endl;root=Euler_method(t);print(sum_num);cout<<"Here is the root of Multistep methods:"<<setw(20)<<setprecision(8)<<root<<endl;cout<<"******************************************************************"<<endl<<endl;}return 0;
}

数值计算一阶常微分方程求解实现相关推荐

2.0.高等数学3-一阶常微分方程求解
一阶常微分方程求解问题引入一阶微分方程的几种情形可分离变量方程齐次微分方程解齐次微分方程例3 例4 一阶线性微分方程标准形式一阶非齐次线性微分方程通解公式伯努利方程三级目录问题引入 ...
matlab欧拉法截断误差,一阶常微分方程欧拉法与梯形公式局部截断误差与p阶精度Range.PPT...
一阶常微分方程欧拉法与梯形公式局部截断误差与p阶精度Range * 一阶常微分方程欧拉法与梯形公式局部截断误差与p阶精度 Range-Kutta公式常微分方程MATLAB求解 <数值分析& ...
一阶欧拉近似matlab,MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程.doc
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程姓名:樊元君学号:2012200902 日期:2012.11.06 一.实验目的掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题 ...
python求解一阶常微分方程
一.用python求解一阶常微分方程 1.求解微分方程需要用到scipy库,pycharm中安装即可,同时需要导入numpy库和matplotlib两个库 2.使用scipy.integrate.od ...
经典龙格-库塔法（四阶龙格-库塔法）求解求一阶常微分方程相应的特解的Python程序
基本原理例题代码 #四阶龙格-库塔法 #求一阶常微分方程,相应的特解 #x变量的区间 a = 0 b = 1 #已知条件 X = [0] Y = [1] h = 0.2 #设置步长 n = (b- ...
一阶常微分方程的数值解法（二阶显式、隐式 Adams 公式及 Milne 方法）
一阶常微分方程的数值解法这里我们介绍二阶显式.隐式 Adams 公式及 Milne 方法求解方程. 题目: 对初值问题 u ′ = u − t 2 , 0 ≤ t ≤ 1 , u ( 0 ) = 0 ...
一阶常微分方程（一）| 存在与唯一性定理 + 变量分离 + 齐次方程
一.一阶常微分方程解的的存在与唯一性定理导数已解出的一阶常微分方程可以表示为如下的一般形式: {dydx=f(x,y)y(x0)=y0(1)\begin{cases} \frac{dy}{dx}=f ...
实验八一阶常微分方程初值问题Matlab实现
不实名感谢某位大兄弟的贡献 [作业内容]: 用Euler法.改进的欧拉法.四阶龙格-库塔法求如下一阶常微分方程 {y′=y−2xyy(0)=1\left\{\begin{array}{l} {y}^{ ...
计算机科学泰勒级数,一阶常微分方程泰勒级数解法的计算机实现.pdf
一阶常微分方程泰勒级数解法的计算机实现第2l卷第西安冶金建筑学院学报 V01．2lNo．4 1989年 12月 J．Xian Inst．ofMetal1．& Cons． ...

数值计算一阶常微分方程求解实现

数值计算一阶常微分方程求解实现相关推荐

最新文章

热门文章