euler.c

用Euler法计算常微分方程数值解

#include<stdio.h>
#include<math.h>
#define N 20

main()
{
	float h;//存放步长
	float a,b;//存放区间端点值
	int n;
	float t0,u0;//存放初值
	float u[N];//存放结果

	void input(float *a,float *b,float *t0,float *u0,float *h);//从文件输入a,b,t0,u0,h
	int get_n(float a,float b,float h);//得到n
	void Euler(float t0,float u0,float h,int n,float u[]);//返回结果，并存放于r[N]
	void output(float t0,float h,float u[],int n);//输出结果

	input(&a,&b,&t0,&u0,&h);

	printf("%f,%f,%f,%f,%f",a,b,t0,u0,h);//测试输入
	printf("\n");

	//n=get_n(a,b,h);
	n=10;
	Euler(t0,u0,h,n,u);
	output(t0,h,u,n);
}

//从文件输入a,b,t0,u0,h
void input(float *pa,float *pb,float *pt0,float *pu0,float *ph)
{
	FILE *f_a_b,*f_t0_u0,*f_h;
	f_a_b=fopen("a_b.dat","r");
	f_t0_u0=fopen("t0_u0.dat","r");
	f_h=fopen("h.dat","r");

	fscanf(f_a_b,"%f %f",pa,pb);
	fscanf(f_t0_u0,"%f %f",pt0,pu0);
	fscanf(f_h,"%f",ph);

	fclose(f_a_b);
	fclose(f_t0_u0);
	fclose(f_h);

}

//计算结果，以数组的形式返回
void Euler(float t0,float u0,float h,int n,float u[])
{
	int i;

	float f(float tn,float un);//计算f(tn,un)
	
	u[0]=u0;
	for(i=1;i<=n;i++)
	{
		u[i]=(4*u[i-1]+h*f(0,u[i-1]))/(4-h);
	}
	printf("\n");
}

//计算f(tn,un)
float f(float tn,float un)
{
	printf("%.2f	",1.0/2.0*un);//测试计算结果
	return 1.0/2.0*un;
}

//输出结果
void output(float t0,float h,float u[],int n)
{
	int i;

	for(i=0;i<=n;i++)
	{
		printf("%.2f  ",t0+i*h);
	}
	printf("\n");
	for(i=0;i<=n;i++)
	{
		printf("%.2f  ",u[i]);
	}
	printf("\n");
	for(i=0;i<=n;i++)
	{
		printf("%.2f  ",exp((t0+i*h)/2));
	}
	printf("\n");
}