数值积分.cpp

来自「用来计算航天器再入质点弹道」· C++ 代码 · 共 82 行

#include"Model.h"
#include "AirModel.h"
//********************************************************************************************
void Euler(int n,double h,double* x,double* f)
{
    Model(x,f);
	for(int	i=0;i<n;i++)
		x[i]+=h*f[i];
}
//**********************************************************************************************
void RK2(int n,double h,double *x,double* f)
{
	double*		xx	=	new double[n];
	double*		k	=	new double[n];
	
	for(int i=0;i<n;i++)
		xx[i]=x[i];
	
	Model(xx,f);
	for(i=0;i<n;i++)
		k[i]=f[i];	    		
	
	for(i=0;i<n;i++)
		xx[i]+=h*k[i];
	
	Model(xx,f);
	for(i=0;i<n;i++)			
	{	
		k[i]+=f[i];     
		x[i]+=h*k[i]/2;
	}
	
	
	delete [] k;
	delete [] xx;
}

//***************************************************************************************************
void RK4(int n,double h,double *x,double* f)
{
	
	double*		k1	=	new double[n];
	double*		k2	=	new double[n];
	double*		k3	=	new double[n];
	double*		k4	=	new double[n];
	double*		xx	=	new double[n];
	
    for(int i=0;i<n;i++)
		xx[i]=x[i];
	
	Model(xx,f);
	for(i=0;i<n;i++)
		k1[i]=f[i];            
	
	for(i=0;i<n;i++)			 
		xx[i]=x[i]+h*k1[i]/2;
	Model(xx,f);
	for(i=0;i<n;i++)					
		k2[i]=f[i];           
	
	for(i=0;i<n;i++)			 
		xx[i]=x[i]+h*k2[i]/2;
	Model(xx,f);
	for(i=0;i<n;i++)					
		k3[i]=f[i];             
	
	for(i=0;i<n;i++)			 
		xx[i]=x[i]+h*k3[i];
	Model(xx,f);
	for(i=0;i<n;i++)					
	{
		k4[i]=f[i];            
		x[i]+=h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6;
	}
	
	delete [] k1;
	delete [] k2;
	delete [] k3;
	delete [] k4;
	delete [] xx;
}