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

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

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

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

赶快添加!!我正好要用到啊
页: [1]
查看完整版本: 一阶常微分方程数初值问题的数值解法(欧拉法)的C++实现