当前位置: 代码迷 >> C语言 >> 我做了一个龙格库塔程序,但不知道对不对!
  详细解决方案

我做了一个龙格库塔程序,但不知道对不对!

热度:185   发布时间:2006-04-27 08:13:00.0
我做了一个龙格库塔程序,但不知道对不对!

#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]);
}


搜索更多相关的解决方案: 龙格库塔  

----------------解决方案--------------------------------------------------------
大哥 ,自己弄个方程运行看看不就直到拉吗??
----------------解决方案--------------------------------------------------------
如果真的想和别人讨论
就加注释,然后说说你的编程思想
----------------解决方案--------------------------------------------------------
如果自信一点呢,就写上[原创]××××××××
如果自谦一点呢,就写上[交流]××××××××
如果自卑一点呢,就写上[讨论]××××××××
但是讨论或求助应顺从3楼提出的合理要求
----------------解决方案--------------------------------------------------------
提示: 作者被禁止或删除 内容自动屏蔽

2006-06-14 21:31:33
cordier

等 级:论坛游民
威 望:1
帖 子:449
专家分:14
注 册:2006-2-9
  得分:0 

怎么说?


----------------解决方案--------------------------------------------------------
提示: 作者被禁止或删除 内容自动屏蔽

2006-06-15 02:48:58
cordier

等 级:论坛游民
威 望:1
帖 子:449
专家分:14
注 册:2006-2-9
  得分:0 
seta[0]我想将来将seta[1]的弧度算为角度输出。因为我们平时都是用角度来看。
给我们一个弧度,我们根本无法了解它是多少。比如弧度0.1,你会觉得它是5度那么大吗?

----------------解决方案--------------------------------------------------------
提示: 作者被禁止或删除 内容自动屏蔽

2006-06-15 03:03:48
cordier

等 级:论坛游民
威 望:1
帖 子:449
专家分:14
注 册:2006-2-9
  得分:0 
不好意思
----------------解决方案--------------------------------------------------------
  相关解决方案