【工程数学】若干种解常微分方程的算法

// ConsoleAppDifferential_Euler_Solu.cpp : 定义控制台应用程序的入口点。
//
/*
*函数功能：欧拉法解常微分方程
*函数原形：void Differential_Euler_Solu(double a,double b)
*参数：double a,double b：所要计算的区间
*返回值：无
*时间复杂度：O(n)
*备注：此法解微分方程精度较低
*日期：2014/12/12
*原创：是
*作者：EbowTang
*Email：tangyibiao520@163.com
*/

#include "stdafx.h"
#include "iostream"
#include <iomanip>

using namespace std;

double Func(double x,double y);
double RealFunc(double x,double y) ;
void Differential_Euler_Solu(double ,double);

int _tmain(int argc, _TCHAR* argv[])
{
cout<<"此微分方程为："<<endl;
cout<<"式子一：y'=y/x-2*y*y"<<endl;
cout<<"式子二：y(0)=0"<<endl;
Differential_Euler_Solu(0,3);
system("pause");
return 0;
}

double Func(double x,double y)
{  if (x!=0.0)
return (y/x-2*y*y);
else
return 1;
}

double RealFunc(double x)//此式子为该微分方程的真实解形式
{
return x/(1+x*x);
}

void Differential_Euler_Solu(double a,double b)
{
double oldX=0.0,oldY=0.0;
double newX=0.0,newY=0.0;

double h=0.2;//步长
for ( int k = 0 ; k <=(b-a)/h ; k++ )
{
oldY = newY;
oldX = newX;
newY = oldY+h*Func(oldX,oldY);//欧拉法的迭代形式
newX = oldX+h;//x递增

cout<<"........................"<<endl;
cout<<"第"<<k+1<<"次迭代计算的结果为"<<endl;
cout<<"x="<<setprecision(5)<<oldX;
cout<<"  真实y="<<setprecision(5)<<RealFunc(oldX);
cout<<"  近似y="<<setprecision(5)<<oldY;
cout<<"  此次误差为="<<setprecision(5)<<oldY-RealFunc(oldX)<<endl;
}
}

<pre name="code" class="cpp">// ConsoleAppDifferential_Improve_Euler_Solu.cpp : 定义控制台应用程序的入口点。
//
/*
*函数功能：改进的欧拉法解常微分方程
*函数原形：void Differential_Improve_Euler_Solu(double a,double b)
*参数：double a,double b：所要计算的区间
*返回值：无
*时间复杂度：O(n)
*备注：此法解微分方程精度相比欧拉法精度明显提高
*日期：2014/12/12
*原创：是
*作者：EbowTang
*Email：tangyibiao520@163.com
*/

#include "stdafx.h"
#include "iostream"
#include <iomanip>

using namespace std;

double Func(double x,double y);
double RealFunc(double x,double y) ;
void Differential_Improve_Euler_Solu(double ,double);

int _tmain(int argc, _TCHAR* argv[])
{
cout<<"此微分方程为："<<endl;
cout<<"式子一：y'=y/x-2*y*y"<<endl;
cout<<"式子二：y(0)=0"<<endl;
Differential_Trapezoid_Solu(0,3);
system("pause");
return 0;
}

double Func(double x,double y)
{  if (x!=0.0)
return (y/x-2*y*y);
else
return 1;
}

double RealFunc(double x)//此式子为该微分方程的真实解形式
{
return x/(1+x*x);
}

void Differential_Improve_Euler_Solu(double a,double b)
{
double oldX=0.0,oldY=0.0;
double newX=0.0,newY=0.0;

double h=0.2;//步长
for ( int k = 0 ; k <=(b-a)/h ; k++ )
{
oldY = newY;
oldX = newX;
newY = oldY+h*Func(oldX,oldY);
newX = oldX+h;//x递增
newY = oldY+(h/2)*(Func(oldX,oldY)+Func(newX,newY));//梯形法的迭代形式

cout<<"........................"<<endl;
cout<<"第"<<k+1<<"次迭代计算的结果为"<<endl;
cout<<"x="<<setprecision(5)<<oldX;
cout<<"  真实y="<<setprecision(5)<<RealFunc(oldX);
cout<<"  近似y="<<setprecision(5)<<oldY;
cout<<"  此次误差为="<<setprecision(5)<<oldY-RealFunc(oldX)<<endl;
}
}

<pre name="code" class="cpp">// ConsoleAppRungeKuttaSolu.cpp : 定义控制台应用程序的入口点。
//
/*
*函数功能：三阶以及四阶RungeKutta法解常微分方程
*函数原形：void Differential_RungeKutta_Solu(double a,double b)
*参数：double a,double b：所要计算的区间
*返回值：无
*时间复杂度：O(n)
*备注：此法解微分方程精度相比欧拉法精度提高
*日期：2014/12/12
*原创：是
*作者：EbowTang
*Email：tangyibiao520@163.com
*/

#include "stdafx.h"
#include "iostream"
#include <iomanip>

using namespace std;

double Func(double x,double y);
double RealFunc(double x,double y) ;
void Differential_RungeKutta_Solu(double ,double);

int _tmain(int argc, _TCHAR* argv[])
{
cout<<"此微分方程为："<<endl;
cout<<"式子一：y'=y/x-2*y*y"<<endl;
cout<<"式子二：y(0)=0"<<endl;
cout<<"/*********************四阶经典龙格法的华丽分割线*******************/"<<endl;
Differential_RungeKutta_Solu(0,3);
system("pause");
return 0;
}

double Func(double x,double y)
{  if (x!=0.0)
return (y/x-2*y*y);
else
return 1;
}

double RealFunc(double x)//此式子为该微分方程的真实解形式
{
return x/(1+x*x);
}

void Differential_RungeKutta_Solu(double a,double b)
{
double oldX=0.0,oldY=0.0;
double newX=0.0,newY=0.0;

double K1=0.0,K2=0.0,K3=0.0,K4;
double h=0.2;//步长
for ( int k = 0 ; k <=(b-a)/h ; k++ )
{
oldY = newY;
oldX = newX;
/************经典四阶龙格法迭代形式*********************/
K1=Func(oldX,oldY);
K2=Func(oldX+h/2,oldY+h/2*K1);
K3=Func(oldX+h/2,oldY+h/2*K2);
K4=Func(oldX+h,oldY+h*K3);
newY=oldY+h/6*(K1+2*K2+2*K3+K4);
newX = oldX+h;//x递增

/******经典三阶龙格法迭代形式
K1=Func(oldX,oldY);
K2=Func(oldX+h/2,oldY+h/2*K1);
K3=Func(oldX+h,oldY-h*K1+2*h*K2);
newY=oldY+h/6*(K1+4*K2+K3);
newX = oldX+h;//x递增
*********************************/
cout<<"第"<<k+1<<"次迭代计算的结果为"<<endl;
cout<<"x="<<setprecision(5)<<oldX;
cout<<"  真实y="<<setprecision(5)<<RealFunc(oldX);
cout<<"  近似y="<<setprecision(5)<<oldY;
cout<<"  此次误差为="<<setprecision(5)<<oldY-RealFunc(oldX)<<endl;
}
}