hjyjxz.c.c

用c语言实现的龙格库塔算法和阿当姆斯算法。是学习计算方法的入门级程序实例。

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

double f;
int N=30;
double h=2,a=0,b=40;
double y[31],x[31],p[31],c[31],m[31],z[31],acc[31];

double hanshu(double x,double z)
{
  double f;
  f=(3150-0.039*z*z)*9.81/(1350-18*x)-9.81;
   return f;
}

void rungerkutta(void)
{
  int i;
  double k1,k2,k3,k4;
  double t=0,m=z[0]=0;
  acc[0]=hanshu(0,0);
  for(i=1;i<=3;i++)
     {
	k1=hanshu(t,m);

	m=z[i-1]+h*k1/2;
	k2=hanshu(t+h/2,m);

	m=z[i-1]+h*k2/2;
	k3=hanshu(t+h/2,m);

	m=z[i-1]+h*k3;
	k4=hanshu(t+h,m);
	y[i]=y[i-1]+h*z[i-1]+h*h*(k1+k2+k3)/6;
	z[i]=z[i-1]+h*(k1+2*k2+2*k3+k4)/6;

	t=i*h;
	acc[i]=hanshu(t,z[i]);
     }
}

void adams(void)
{
  int i,j;
  double k[5],t;
  double r[4];
  r[0]=z[0];r[1]=z[1];r[2]=z[2];r[3]=z[3];

   for(j=0;j<4;j++)
   {
    t=j*h;
    k[j]=hanshu(t,r[j]);
   }
   for(i=4;i<=N;i++)
     {
	t=t+h;
	p[i]=r[3]+h*(55*k[3]-59*k[2]+37*k[1]-9*k[0])/24;
	m[i]=p[i]+(251.0/270)*(c[i-1]-p[i-1]);
	k[4]=hanshu(t,m[i]);
	c[i]=r[3]+h*(9*k[4]+19*k[3]-5*k[2]+k[1])/24;
	z[i]=c[i]-19.0/270*(c[i]-p[i]);
	k[4]=hanshu(t,z[i]);
	r[3]=z[i];

	acc[i]=hanshu(t,z[i]);
	y[i]=y[i-2]+h/3.0*(z[i-2]+4.0*z[i-1]+z[i]);
	for(j=0;j<4;j++)
	  k[j]=k[j+1];
     }
}

void output(void)
{
  int i;
  for(i=0;i<=N;i++)
   {
    x[i]=h*i;

   }

  printf("k     x[k]         y[k]       v[k]        a[k]           \n");
  for(i=0;i<=3;i++)
    printf("%-5d%-.4f      %3.6e      %3.6e     %.5e    \n",i,x[i],y[i],z[i],acc[i]);
  for(i=4;i<=N;i++)
    printf("%-5d%-.4f     %3.6e    %3.6e     %.5e \n",i,x[i],y[i],z[i],acc[i]);
}

void main()
{
   rungerkutta();
   adams();
   output();
}