程序一直运行但不出结果，请教原因。

bingxue

#include<iostream>
#include<fstream>
#include<cmath>
using namespace std;
ofstream rotor;
double w=874.8*2*3.14/60,e=0.0003;
double m=420,D=0.114,L=0.0285,R=D/2,T=R/2;
double k=2.104e8,c=9300,ka=k/10,kb=k/15;
double pai=3.14,g=9.8;
double A1=e*w*w;
double a=acos((R-T)/R);
double b(double t,double p);
double fb(double t,double p);
double kx(double t,double p);
double kxy(double t,double p);
double ky(double t,double p);
double h1(double X,double x,double y,double t,double p);
double h2(double Y,double x,double y,double t,double p);
void main()
{
  double X,Y,x,y;
  double K1,L1,M1,N1;
  double K2,L2,M2,N2;
  double K3,L3,M3,N3;
  double K4,L4,M4,N4;
  rotor.open("e:\\rotor.txt");
  double h=0.001,t=0;
  X=Y=0,x=y=0;
  for(int i=1;i<30000;i++)
{
   t=t+h;
   double p=atan(y/x);
   K1=X;
   L1=h1(X,x,y,t,p);
   M1=Y;
   N1=h2(Y,x,y,t,p);

   p=atan((y+h*M1/2)/(x+h*K1/2));
   K2=X+h*L1/2;
   L2=h1(X+h*L1/2,x+h*K1/2,y+h*M1/2,t+h/2,p);
   M2=Y+h*N1/2;
   N2=h2(Y+h*N1/2,x+h*K1/2,y+h*M1/2,t+h/2,p);

   p=atan((y+h*M2/2)/(x+h*K2/2));  
   K3=X+h*L2/2;
   L3=h1(X+h*L2/2,x+h*K2/2,y+h*M2/2,t+h/2,p);
   M3=Y+h*N2/2;
   N3=h2(Y+h*N2/2,x+h*K2/2,y+h*M2/2,t+h/2,p);

   p=atan((y+h*M3)/(x+h*K3));
   K4=X+h*L3;
   L4=h1(X+h*L3,x+h*K3,y+h*M3,t+h,p);
   M4=Y+h*N3;
   N4=h2(Y+h*N1,x+h*K1,y+h*M1,t+h,p);

   x=x+h*(K1+2*K2+2*K3+K4)/6;
   X=X+h*(L1+2*L2+2*L3+L4)/6;
   y=y+h*(M1+2*M2+2*M3+M4)/6;
   Y=Y+h*(N1+2*N2+2*N3+N4)/6;
   rotor<<x<<endl;
}
rotor.close();
}
double b(double t,double p)
{
  double n=w*t-p+pai/2;
  if(n<2*pai)  return n;
  else
     for (int j=1;n>=2*pai;j++)
    { n=n-2*j*pai;}
    return n;
}
double fb(double t,double p)
{
   if(b(t,p)<=(pai/2-a)||b(t,p)>=(-pai/2+a))  return 1;
   if(b(t,p)<=(pai/2+a)||b(t,p)>=(pai/2-a))  return (1+cos((fb(t,p)-pai/2+a)*pai/(2*a)))/2;
   if(b(t,p)<=(3*pai/2-a)||b(t,p)>=(pai/2+a))  return 0;
   if(b(t,p)<=(3*pai/2+a)||b(t,p)>=(3*pai/2-a))  return (1+cos((fb(t,p)-3*pai/2+a)*pai/(2*a)))/2;
}
double kx(double t,double p)
{
  return k-fb(t,p)*(ka*cos(w*t)*cos(w*t)+kb*sin(w*t)*sin(w*t));
}
double kxy(double t,double p)
{
  return -fb(t,p)*(ka-kb)*cos(w*t)*cos(w*t)*sin(w*t)*sin(w*t);
}
double ky(double t,double p)
{
  return k-fb(t,p)*(ka*sin(w*t)*sin(w*t)+kb*cos(w*t)*cos(w*t));
}
double h1(double X,double x,double y,double t,double p)
{
  return -c*X/m-kx(t,p)*x/m-kxy(t,p)*y/m+A1*cos(w*t);
}
double h2(double Y,double x,double y,double t,double p)
{
  return -g-c*Y/m-ky(t,p)*y/m-kxy(t,p)*x/m+A1*sin(w*t);
}


谢谢！

kangarooli

不知道你程序有什么问题，我也遇到过，我的感觉是数值太复杂，所以造成运算不出来

tokyxp

我不知道这段程序用途是什么，但是这段 程序编译没有问题，
但程序中存在问题
1、 Math.Atan((y + h * M1 / 2) / (x + h * K1 / 2)); 分母为0 程序返回值为nan，不能够计算。
2、程序中有返回值语句，部分返回值没有设定
        double fb(double t, double p)
        {
            if (b(t, p) <= (pai / 2 - a) || b(t, p) >= (-pai / 2 + a)) return 1;
            if (b(t, p) <= (pai / 2 + a) || b(t, p) >= (pai / 2 - a)) return (1 + Cos((fb(t, p) - pai / 2 + a) * pai / (2 * a))) / 2;
            if (b(t, p) <= (3 * pai / 2 - a) || b(t, p) >= (pai / 2 + a)) return 0;
            if (b(t, p) <= (3 * pai / 2 + a) || b(t, p) >= (3 * pai / 2 - a)) return (1 + Cos((fb(t, p) - 3 * pai / 2 + a) * pai / (2 * a))) / 2;
        }
如果所有的if条件都不满足的话，你说它能够返回什么呢？