ode.cpp

基于龙格库塔算法的矩阵微分方程组求解子程序

///////本程序是在求解电路暂态分析时想到的，可以求解多阶的微分方程组，特别是基于矩阵形式的微分方程组
void get_ode(float **a,float*m)
{
	float h;
	float	*k1,*k2,*k3,*k4;
	extern int lc;
	extern float **y,ti;
	int i,j,k;
	h=ti/1000;
	k1=new float[lc];
	k2=new float[lc];
	k3=new float[lc];
	k4=new float[lc];
    for(i=0;i<lc;i++)
	{  
		y[i][0]=state[i];
	}
	for(k=0;k<999;k++)
{
		for(i=0;i<lc;i++)
		{
		k1[i]=0;
	    k2[i]=0;
		k3[i]=0;
		k4[i]=0;
		}
	for(i=0;i<lc;i++)
	{
		for(j=0;j<lc;j++)
	      k1[i]+=a[i][j]*y[j][k];
		k1[i]+=m[i];
	}////////////计算k1值在x处导数
		for(i=0;i<lc;i++)
		{
		for(j=0;j<lc;j++)
			k2[i]+=a[i][j]*(k1[j]*h/3+y[j][k]);
		k2[i]+=m[i];
		}////////计算k2值在x+h/2处的导数
		for(i=0;i<lc;i++)
		{
	    	for(j=0;j<lc;j++)
				k3[i]+=a[i][j]*(k2[j]*h/3+y[j][k]);
            k3[i]+=m[i];
		}/////计算k3值计算在x+h/2处的导数
        	for(i=0;i<lc;i++)
		{
	    	for(j=0;j<lc;j++)
				k4[i]+=a[i][j]*(k3[j]*h+y[j][k]);
            k4[i]+=m[i];
		}
		for(i=0;i<lc;i++)
		y[i][k+1]=y[i][k]+(k1[i]+3*k2[i]+3*k3[i]+k4[i])*h/8;
}
	delete []k1;
		delete []k2;
			delete []k3;
				delete []k4;

}