adams.c

我自己写的数值分析课程设计的源代码

#include<stdio.h>
#include<math.h>

#define n 20//19 20//17 50

double yi =0;//17 1//19 0
double ai =0;
double bi =1;

double (*True) (double) = NULL;
double ERR=0.0;

double x[n+1]={0},p[n+1]={0},m[n+1]={0},c[n+1]={0},y[n+1]={0};
double h,k1,k2,k3,k4;

void	kutta(int m);
void	print(double *a,int m,char *c);
double	f(int m);
double	ff(double a,double b);
double	g(double x);
double	k(double x);


void main()
{
	int i,j;
	FILE *fout;
	h=(bi-ai)/n;
	//*True=g;//17
	True=&k;//19

	for(i=0;i<=n;i++)
		x[i]+=h*i;
	y[0]=yi;
	for(i=1;i<=n;i++)
		kutta(i);
	
	if(fout=fopen("adams.txt","w"))
	{
		fprintf(fout,"Runge-Kutta:\n");
		for(j=0;j<=n;j++)		
		{
			ERR = fabs(y[j]-True(x[j]));
			fprintf(fout,"(%4.2f,%f)    True = %-12g    ERR = %-12g\n",x[j],y[j],True(x[j]),ERR);
		}
		fprintf(fout,"\n");
	}
	else
		printf("can't creat data!\n");
	
	for(i=3;i<=n-1;i++)
	{
		p[i+1]=y[i]+h*(55*f(i)-59*f(i-1)+37*f(i-2)-9*f(i-3))/24;
		m[i+1]=p[i+1]+251*(c[i]-p[i])/270;
		c[i+1]=y[i]+h*(9*ff(x[i+1],m[i+1])+19*f(i)-5*f(i-1)+f(i-2))/24;
		y[i+1]=c[i+1]-19*(c[i+1]-p[i+1])/270;
	}
	
	fprintf(fout,"Adams:\n");
	for(j=0;j<=n;j++)		
	{
		ERR = fabs(y[j]-True(x[j]));
		fprintf(fout,"(%4.2f,%f)    True = %-12g    ERR = %-12g\n",x[j],y[j],True(x[j]),ERR);
	}
	fprintf(fout,"\n");
	fclose(fout);
}
	


void print(double *a,int m,char *c)
{
	int i;
	FILE *fout;
	if(fout=fopen("adams.txt","a+"))
	{
		for(i=0;i<=m;i++)
		{
			fprintf(fout,"%s[%d]=%f\n",c,i,a[i]);
		}
		fprintf(fout,"\n");
		fclose(fout);
	}
	else
		printf("Disk is full!\n");
}

double f(int m)
{
	//return x[m]+y[m];//17
	return 2*pow(x[m],3)-2*x[m]*y[m];//19
}

double ff(double a,double b)
{
	//return a+b;//17
	return 2*pow(a,3)-2*a*b;//19
}

void kutta(int m)
{
	k1=h*ff(x[m-1],y[m-1]);
	k2=h*ff(x[m-1]+.5*h,y[m-1]+.5*k1);
	k3=h*ff(x[m-1]+.5*h,y[m-1]+.5*k2); 
	k4=h*ff(x[m-1]+h,y[m-1]+k3);
	y[m]=y[m-1]+(k1+k2*2+k3*2+k4)/6;
	
}

double g(double x)//17
{
	return 2*exp(x)-x-1;
}

double k(double x)//19
{
	return exp(-pow(x,2))+pow(x,2)-1;
}