/*以下消除了楼主源程序编译时的5个warning*/
/*至于其专业功能正常与否需要楼主自行判断*/
#include <math.h>
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
main()
{
long i=1,n=2;//原先楼主错将n定义为double型
double omga,l,g=9.8,t=0,beta,F,miu; //A,deta冗余
double seta[3],dydx[3],yout[3],h,dt=0.01;//n移出
FILE *fp; //fpdata冗余
void function(double t,double seta[],double dydx[],double omga,double beta,double F,double miu);
void rk4(double y[],double yn[],int n,double t,double h,double yout[],double omga,double beta,double F,double miu);
system("cls");
l=1.0;
omga=sqrt(g/l);
beta=0.25;
miu=2.0/3;
seta[1]=1;
seta[2]=0;
F=1.15;

if((fp=fopen("seta.txt","w"))==NULL)
{
printf("the file:seta.txt can not open:Press any key to exit:");
getch();
exit(1);
}
function(t,seta,dydx,omga,beta,F,miu);
h=dt;
fprintf(fp,"%.5lf\t%.5lf\t%.5lf\t%.5lf\n",t,seta[1],seta[2],miu*t);
do
{
function(t,seta,dydx,omga,beta,F,miu);
t+=dt;
rk4(seta,dydx,n,t,h,yout,omga,beta,F,miu);
seta[1]=yout[1];
seta[2]=yout[2];
i++;
fprintf(fp,"%.5lf\t%.5lf\t%.5lf\t%.5lf\n",t,yout[1],yout[2],miu*t);
}
while(!kbhit());
getch();
printf("F=%.5lf\tend\nPress any key to continue:",F);
fclose(fp);
printf("i=%ld\n",i);
printf("beta=%.7lf\nmiu=%.7lf\nF=%.7lf\nomga=%.7lf\n",beta,miu,F,omga);
getch();
}

void rk4(double y[],double yn[],int n,double t,double h,double yout[],double omga,double beta,double F,double miu)
{
int i;
double ytemp[3],k2[3],k3[3],k4[3],th,h2,h6;
h2=h*0.5;
h6=h*1.0/6;
th=t+h2;
for(i=1;i<3;i++)
{
ytemp[i]=0;
k2[i]=0;
k3[i]=0;
k4[i]=0;
}
for(i=1;i<=n;i++)
ytemp[i]=y[i]+h2*yn[i];
function(th,ytemp,k2,omga,beta,F,miu);
for(i=1;i<=n;i++)
ytemp[i]=y[i]+h2*k2[i];
function(th,ytemp,k3,omga,beta,F,miu);
for(i=1;i<=n;i++)
ytemp[i]=y[i]+h*k3[i];
function(t+h,ytemp,k4,omga,beta,F,miu);
for(i=1;i<=n;i++)
yout[i]=y[i]+h6*(yn[i]+2.0*k2[i]+2.0*k3[i]+k4[i]);
}

void function(double t,double seta[],double dydx[],double omga,double beta,double F,double miu)
{
dydx[1]=seta[2];
dydx[2]=-omga*omga*sin(seta[1]);
}

[此贴子已经被作者于2006-6-16 18:28:37编辑过]

[转载]供楼主参考
//分别用2、3、4阶龙格库塔方法
//求解y'=y-2*x/y(0<x<1); y0=1;

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

typedef double R8;
void RungeKutta(R8,R8,R8,int,R8*,R8*,int,R8(*function)(R8,R8));
R8 function(R8,R8);

main()
{ int i;
R8 x[6],y[6];
printf("用二阶龙格-库塔方法\n");
RungeKutta(1,0,1,5,x,y,2,function);
for(i=0;i<6;i++)
printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);
printf("用三阶龙格-库塔方法\n");
RungeKutta(1,0,1,5,x,y,3,function);
for(i=0;i<6;i++)
printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);
printf("用四阶龙格-库塔方法\n");
RungeKutta(1,0,1,5,x,y,4,function);
for(i=0;i<6;i++)
printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);
}

void RungeKutta(R8 y0,R8 a,R8 b,int n,R8 *x,R8 *y,int style,R8(*function)(R8,R8))
{ /* n表示等分数,n+1表示头尾总数 */
R8 h=(b-a)/n,k1,k2,k3,k4;
int i;
x[0]=a;
y[0]=y0;
switch(style)
{
case 2:
for(i=0;i<n;i++)
{
x[i+1]=x[i]+h;
k1=function(x[i],y[i]);
k2=function(x[i]+h/2,y[i]+h*k1/2);
y[i+1]=y[i]+h*k2;
}
break;

case 3:
for(i=0;i<n;i++)
{
x[i+1]=x[i]+h;
k1=function(x[i],y[i]);
k2=function(x[i]+h/2,y[i]+h*k1/2);
k3=function(x[i]+h,y[i]-h*k1+2*h*k2);
y[i+1]=y[i]+h*(k1+4*k2+k3)/6;
}
break;

case 4:
for(i=0;i<n;i++)
{
x[i+1]=x[i]+h;
k1=function(x[i],y[i]);
k2=function(x[i]+h/2,y[i]+h*k1/2);
k3=function(x[i]+h/2,y[i]+h*k2/2);
k4=function(x[i]+h,y[i]+h*k3);
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;
}
break;
}
}

R8 function(R8 x,R8 y)//实际是导函数
{
return y-2*x/y;
}

花了好长时间，终于得到楼主程序的简化版
/*++++++++++++++++++++++++++++++++++++++++
楼主所研究的是：摆长＝1米的单摆的运动轨迹。
t＝0时单摆偏离平衡位置1弧度(57.3度)初速＝0
龙格–库塔运行结果周期Ｔ＝2.14秒(取g＝9.8)
该结果与单摆小幅摆动时的周期2.00秒略有不同
+++++++++++++++++++++++++++++++++++++++++*/
#include<math.h>
#include<conio.h>
#include<stdio.h>
#include<stdlib.h>

double omga;

void RK4(double *y,double *yn,double t,double h,double *yout);

main()
{
long i=0;
char key;
double seta[2]={1,0},dydx[2],yout[2],t=0,dt=0.01;
double l=1,g=9.8; //单摆长度l=1米,g=9.8(重力加速度)
omga = sqrt(g/l); //单摆角频率ω
dydx[0]=seta[1];
dydx[1]=-omga*omga*sin(seta[0]);
printf("时刻(sec)\t张角(arc)\t速度(m/s)\n");
printf("%lf\t%lf\t%lf\n",t,seta[0],seta[1]);
while(1)
{
if(++i%24==0)key=getch();
if(key==0x1b)break;
dydx[0]=seta[1];
dydx[1]=-omga*omga*sin(seta[0]);
t+=dt;
RK4(seta,dydx,t,dt,yout);
seta[0]=yout[0];
seta[1]=yout[1];
printf("%lf\t%lf\t%lf\n",t,yout[0],yout[1]);
}
printf("i=%ld\n",i);
printf("omga=%lf\n",omga);
}

void RK4(double *y,double *yn,double t,double h,double *yout)
{
double ytemp[2],k2[2],k3[2],k4[2];
ytemp[0]=y[0]+h/2*yn[0];
ytemp[1]=y[1]+h/2*yn[1];
k2[0]=ytemp[1];
k2[1]=-omga*omga*sin(ytemp[0]);
ytemp[0]=y[0]+h/2*k2[0];
ytemp[1]=y[1]+h/2*k2[1];
k3[0]=ytemp[1];
k3[1]=-omga*omga*sin(ytemp[0]);
ytemp[0]=y[0]+h*k3[0];
ytemp[1]=y[1]+h*k3[1];
k4[0]=ytemp[1];
k4[1]=-omga*omga*sin(ytemp[0]);
yout[0]=y[0]+h/6*(yn[0]+2*k2[0]+2*k3[0]+k4[0]);
yout[1]=y[1]+h/6*(yn[1]+2*k2[1]+2*k3[1]+k4[1]);
}

[此贴子已经被作者于2006-6-18 8:57:30编辑过]