#include <math.h>
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
main()
{
long i=1;
double omga,l,g=9.8,t=0,deta,beta,F,miu,A;
double seta[3],dydx[3],yout[3],n=2,h,dt=0.01;
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();
}
/*====================rk4()==================*/
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=1.0*h/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]);
}
/*===============function()==================*/
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]);
}
