ruk.cpp

四阶龙格库达法求解微分方程组

// Ruk.cpp: implementation of the CRuk class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Ruk.h"
#include "math.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CRuk::CRuk()
{
	n=3;    //方程的个数
	y[0]=0;   //赋初值
    //^^^^^^^^^^^^^^^^^^^
	y[n-1]=0;

}

CRuk::~CRuk()
{

}
void CRuk::dif()          //微分方程组
{
	d[0] = m_dV_T*cos(m_dSeita_T)*cos(m_dFai_T-y[2])*cos(y[1])+m_dV_T*sin(m_dSeita_T)*sin(y[1])
		   -m_dV_M*cos(m_dSeita_M)*cos(m_dFai_M-y[2])+m_dV_M*sin(m_dSeita_M)*sin(y[1]);
	d[1] = (m_dV_T*cos(m_dSeita_T)*cos(m_dFai_T-y[2])*sin(y[1])+m_dV_T*sin(m_dSeita_T)*cos(y[1])
		   +m_dV_M*cos(m_dSeita_M)*cos(m_dFai_M-y[2])*sin(y[1])-m_dV_M*sin(m_dSeita_M)*cos(y[1]))/y[0];
	d[2] = (m_dV_T*cos(m_dSeita_T)*sin(m_dFai_T-y[2])+m_dV_M*cos(m_dSeita_M)*sin(m_dFai_M-y[2]))/(y[0]*cos(y[1]));	
}

void CRuk::ruk4()         //龙格库达法求解
{ 
	int i;
	double rk1[5],rk2[5],rk3[5],ys[5],dt02;
	dt02=0.5*dt;
	dif();
	for(i=0;i<n;i++)
	{
		rk1[i]=d[i]*dt02;
		ys[i]=y[i];
		y[i]=ys[i]+rk1[i];
	}	
    dif();
	for(i=0;i<n;i++)
	{
		rk2[i]=d[i]*dt02;
		y[i]=ys[i]+rk2[i];
	}
	dif();
	for(i=0;i<n;i++)
	{
		rk3[i]=d[i]*dt;
		y[i]=ys[i]+rk3[i];
	}
    dif();
	for(i=0;i<n;i++)
	{
		y[i]=ys[i]+(2*(rk1[i]+rk3[i])+4*rk2[i]+d[i]*dt)/6;
	}
}