算法及编程语言 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

christy 发表于 2005-11-25 19:56

一阶常微分方程数初值问题的数值解法（欧拉法）的C++实现

一阶常微分方程一般可以写成如下一般形式： dy/dx=f(x,y) f(x0)=y0 对于此类初值问题，可以用欧拉法求解，欧拉法主要有三种计算格式： 1.显式格式； 2,隐式格式； 3.梯形格式； 写了一个通用的C++类CCDE(Class of Constant Differential Equation)，把以上三种求解格式都包含进去了。 整个程序的源代码由头文件CDE.h和实现文件CDE.CPP组成； 源代码如下： 1.头文件CDE.h //CDE.h，头文件的全部内容： // CDE.h: interface for the CCDE class. // ////////////////////////////////////////////////////////////////////// #if !defined(AFX_CDE_H__712F05EE_06CA_4113_832E_4A8B6E94FE5B__INCLUDED_) #define AFX_CDE_H__712F05EE_06CA_4113_832E_4A8B6E94FE5B__INCLUDED_ #if _MSC_VER > 1000 #pragma once #endif // _MSC_VER > 1000 #include "stdafx.h" typedef struct { double x; double y; } PT;//解的横坐标和纵坐标所代表的点； class CCDE { public: PT* R; CCDE(); virtual ~CCDE(); public: void DrawData(CDC* pDC); void InitData(double x0,double xe,double y0,UINT N); void Euler_Explicit(void); double (*pf)(double x,double y); private: BOOL m_bInit; double y0;//求解区域的起点初值； double x0;//求解区域的起点的横坐标； double xe;//求解区域的终点的横坐标； UINT N;//求解区域平均等份数； double h;//求解区域的划分步长；h由x0,xe,N确定==>h=(xe-x0)/N; public: void Euler_UnExplicit(double error); void Euler_TX(double error); }; #endif // !defined(AFX_CDE_H__712F05EE_06CA_4113_832E_4A8B6E94FE5B__INCLUDED_) 2.实现文件CDE.cpp; // CDE.cpp: implementation of the CCDE class. // ////////////////////////////////////////////////////////////////////// #include "stdafx.h" #include "Ex1.h" #include "CDE.h" #include ".\cde.h" #include "math.h" #ifdef _DEBUG #undef THIS_FILE static char THIS_FILE[]=__FILE__; #define new DEBUG_NEW #endif ////////////////////////////////////////////////////////////////////// // Construction/Destruction ////////////////////////////////////////////////////////////////////// CCDE::CCDE() { pf=NULL; m_bInit=false; R=NULL; } CCDE::~CCDE() { if(R!=NULL) { delete [] R; R=NULL; } } void CCDE::Euler_Explicit() { if(pf==NULL) { AfxMessageBox("****** 函数未指定! ******\n求解前应先使函数指针pf指向函数f(x,y)"); if(R!=NULL) { delete [] R; R=NULL; exit(1); } } if(!m_bInit) { AfxMessageBox("***已知条件未初始化***\n在求解前先调用InitData()函数进行初始化"); if(R!=NULL) { delete [] R; R=NULL; exit(1); } } double fxy; for(UINT i=1;i<=N;i++) { R.x=x0+h*i; fxy=(*pf)(R.x-h,R.y); R.y=R.y+h*fxy; } R.x=xe; } /*---------------------------------------------------------------- x0:求解区域的起点的横坐标； xe:求解区域的终点的横坐标； y0:给定的初值；(未知函数在起点坐标处的值) N:求解区域的平均等份数； /*-----------------------------------------------------------------*/ void CCDE::InitData(double x0, double xe, double y0, UINT N) { this->x0=x0; this->xe=xe; this->y0 =y0; this->N=N; /********************************************/ h=(xe-x0)/N; m_bInit=true; R=new PT; R.x=x0; R.y=y0; R.x=xe; } void CCDE::DrawData(CDC *pDC) { char ch; CString str; for(UINT i=0;i<=N;i++) { sprintf(ch,"(%0.1f , %0.4f)",R.x,R.y); pDC->TextOut(50+90*i,50+20*i,ch); } } void CCDE::Euler_UnExplicit(double error) { if(pf==NULL) { AfxMessageBox("****** 函数未指定! ******\n求解前应先使函数指针pf指向函数f(x,y)"); if(R!=NULL) { delete [] R; R=NULL; exit(1); } } if(!m_bInit) { AfxMessageBox("***已知条件未初始化***\n在求解前先调用InitData()函数进行初始化"); if(R!=NULL) { delete [] R; R=NULL; exit(1); } } double fxy; double oldR; double f; for(UINT i=1;i<=N;i++) { R.x=x0+h*i; fxy=(*pf)(R.x,R.y); R.y=R.y+h*fxy;//先用显式公式估算初值R.y； oldR=R.y;//保存迭代前的值； f=(*pf)(R.x,R.y);//把上面用隐式公式算出的R.y带入公式迭代； R.y=R.y+h*f;//隐式公式； while(fabs(R.y-oldR)>=error)//比较相邻两次的R.y； { oldR=R.y; f=(*pf)(R.x,R.y); R.y=R.y+h*f; } } R.x=xe;//求解区域的终点的横坐标用已知的值代替，以避免使用R.x=x0+h*i带来的误差； } void CCDE::Euler_TX(double error) { if(pf==NULL) { AfxMessageBox("****** 函数未指定! ******\n求解前应先使函数指针pf指向函数f(x,y)"); if(R!=NULL) { delete [] R; R=NULL; exit(1); } } if(!m_bInit) { AfxMessageBox("***已知条件未初始化***\n在求解前先调用InitData()函数进行初始化"); if(R!=NULL) { delete [] R; R=NULL; exit(1); } } double fxy; double oldR; double f1,f2; for(UINT i=1;i<=N;i++) { R.x=x0+h*i; fxy=(*pf)(R.x,R.y); R.y=R.y+h*fxy; oldR=R.y; f1=(*pf)(R.x,R.y); f2=(*pf)(R.x,R.y); R.y=R.y+h*(f1+f2)/2; while(fabs(R.y-oldR)>=error) { oldR=R.y; f1=(*pf)(R.x,R.y); f2=(*pf)(R.x,R.y); R.y=R.y+h*(f1+f2)/2; } } R.x=xe; } 使用方法：只需要把头文件CDE.h和实现文件CDE.cpp加入到工程中，并在用到的地方把头文件include 进来即可。 使用说明： 1.为了程序的通用性，类中使用了函数指针，求解的时候，需要根据实际情况，自己定义函数f(x,y),然后使函数指针pf指向此函数； 2.调用函数InitData ()进行初始化。 3.调用其中一种格式进行计算； 4.计算结果保存在成员变量R中；R[数组索引].x代表横坐标，R[数组索引].y代表对应的纵坐标。 比如:R.x，R.y代表编号为1的点的横坐标和纵坐标值。 注意：编号从0开始，而且题目给定的初始点的编号为0； 下面是示例程序中的片断： ........................................... pDoc->eq.InitData(0,0.5,2,5); pDoc->eq.pf=func; pDoc->eq.Euler_Explicit(); ......................................... 通过对几个具体的算例，把计算结果与Matlab的计算结果进行了对比，计算结果准确可靠。 局限性： 1.因为使用了AfxMessageBox()函数来给出错误提示信息，所以该类不能用于Console类型的程序，有空将提供用于Console的版本。 有空将继续添加通用的龙格－库塔算法以及两点边值问题的算法。 还打算写一个通用的矩阵类，实现常见的矩阵运算。 还有线性方程组的数值解法的C++实现等。

brettlin 发表于 2006-4-14 14:50

赶快添加!!我正好要用到啊

页: [1]

声振论坛's Archiver

一阶常微分方程数初值问题的数值解法（欧拉法）的C++实现