difeguation.cpp

常用数值算法的C++源码

#include "stdio.h"
#include "math.h"
#define  SIZE 100            //定义存储节点的数组大小
double   y_x[SIZE];          //存放真值的数组   
double   y_0[4];             //存放初值的数组
FILE     *fp=fopen("difEguation.txt","w");

//f(x,y)
double f1(double x,double y)
{
	double f = -pow(x,2)*pow(y,2);
	return f;
}

//精确解函数
double f2(double x)
{
	double f = 3/(1+pow(x,3));
	return f;
}

//RK4
double *RK4(double x[],double y0,double h,int n)
{
	int    i;
	double k1,k2,k3,k4;
	double y[SIZE];
	double err[SIZE];
	double *r = y;

	for(i=0,y[0]=y0;i<n-1;i++)
	{
		k1 = f1(x[i],y[i]);
		k2 = f1(x[i]+0.5*h,y[i]+0.5*h*k1);
		k3 = f1(x[i]+0.5*h,y[i]+0.5*h*k2);
		k4 = f1(x[i+1],y[i]+h*k3);
		y[i+1] = y[i]+h*(k1+2*k2+2*k3+k4)/6;
    }
	for(i=0;i<4;i++) y_0[i] = y[i];
    fprintf(fp,"\n    RK4计算结果          误差\n");
	for(i=0;i<n;i++) 
	{
		err[i] = fabs(y[i]-y_x[i]);
		fprintf(fp,"     %f          %f\n",y[i],err[i]);
	}
	return r;
}

//AB4
double *AB4(double x[],double y0,double y1,double y2,double y3,double h,int n)
{
	int    i;
	double y[SIZE];
	double err[SIZE];
	double *r = y;
	y[0] = y0;
	y[1] = y1;
	y[2] = y2;
	y[3] = y3;
	for(i=3;i<n-1;i++)
		y[i+1] = y[i]+h*(55*f1(x[i],y[i])-59*f1(x[i-1],y[i-1])+37*f1(x[i-2],y[i-2])-9*f1(x[i-3],y[i-3]))/24;
	fprintf(fp,"\n    AB4计算结果          误差\n");
	for(i=0;i<n;i++) 
	{
		err[i] = fabs(y[i]-y_x[i]);
		fprintf(fp,"     %f          %f\n",y[i],err[i]);
	}
	return r;
}

//AB4-AM4
double *AB4AM4(double x[],double y0,double y1,double y2,double y3,double h,int n)
{
	int    i;
	double y[SIZE];
	double err[SIZE];
	double *r = y;
	y[0] = y0;
	y[1] = y1;
	y[2] = y2;
	y[3] = y3;
	for(i=3;i<n-1;i++)
	{
		y[i+1] = y[i]+h*(55*f1(x[i],y[i])-59*f1(x[i-1],y[i-1])+37*f1(x[i-2],y[i-2])-9*f1(x[i-3],y[i-3]))/24;
		y[i+1] = y[i]+h*(9*f1(x[i+1],y[i+1])+19*f1(x[i],y[i])-5*f1(x[i-1],y[i-1])+f1(x[i-2],y[i-2]))/24;
	}
	fprintf(fp,"\n  AB4-AM4计算结果        误差\n");
	for(i=0;i<n;i++) 
	{
		err[i] = fabs(y[i]-y_x[i]);
		fprintf(fp,"     %f          %f\n",y[i],err[i]);
	}
	return r;
}

//改进AB4-AM4
double *AB4AM4Impro(double x[],double y0,double y1,double y2,double y3,double h,int n)
{
	int    i;
	double temp;
	double y[SIZE];
	double err[SIZE];
	double *r = y;
	y[0] = y0;
	y[1] = y1;
	y[2] = y2;
	y[3] = y3;
	for(i=3;i<n-1;i++)
	{
		temp   = y[i]+h*(55*f1(x[i],y[i])-59*f1(x[i-1],y[i-1])+37*f1(x[i-2],y[i-2])-9*f1(x[i-3],y[i-3]))/24;
		y[i+1] = y[i]+h*(9*f1(x[i+1],temp)+19*f1(x[i],y[i])-5*f1(x[i-1],y[i-1])+f1(x[i-2],y[i-2]))/24;
		y[i+1] = (251*y[i+1]+19*temp)/270;
	}
	fprintf(fp,"\n改进AB4-AM4计算结果      误差\n");
	for(i=0;i<n;i++) 
	{
		err[i] = fabs(y[i]-y_x[i]);
		fprintf(fp,"     %f          %f\n",y[i],err[i]);
	}
	return r;
}

void main()
{
	double  a = 0;
	double  b = 1.5;
	double  h = 0.1;
    int     n = (int)((b-a)/h+1);          //节点数
	double  x[SIZE];

    for(int i=0;i<n;i++) 
	{
		x[i]   = a+i*h;
		y_x[i] = f2(x[i]);
	}
	RK4(x,3,0.1,n);
	AB4(x,y_0[0],y_0[1],y_0[2],y_0[3],0.1,n);
	AB4AM4(x,y_0[0],y_0[1],y_0[2],y_0[3],0.1,n);
	AB4AM4Impro(x,y_0[0],y_0[1],y_0[2],y_0[3],0.1,n);
	printf("运行结果保存至文件difEguation.txt\n");
}