刚才上传的附件，不方便看。下面是用C编写的runge kutta求解转子振动方程

胡晓宇

不好意思，刚才上传附件，不方便大家看。下面是我用C编的求解轴承转子振动的runge-kutta法，大家来看一下。,为什么我求解了1000个周期，其转子中心位移越来越大呢，是不是结果发散了？请大家多批评，指点。谢谢！
#include "stdlib.h"
#include "math.h"

void rkt2(t,h,y,n,eps,f)
void(*f)();
int n;
double t,h,eps,y[];
{ int m,i,j,k;
  double hh,p,dt,x,tt,q,a[4],*g,*b,*c,*d,*e;
  g=malloc(n*sizeof(double));
  b=malloc(n*sizeof(double));
  c=malloc(n*sizeof(double));
  d=malloc(n*sizeof(double));
  e=malloc(n*sizeof(double));
  hh=h;
  m=1;
  p=1+eps;
  x=t;
  for(i=0; i<=n-1; i++)
      c=y;
  while(p>=eps)
  {
      a[0]=hh/2.0;
      a[1]=a[0];
      a[2]=hh;a[3]=hh;
      for(i=0; i<=n-1;i++)
      {
          g=y;
          y=c;
      }
      dt=h/m;
      t=x;
      for(j=0;j<=m-1; j++)
      {
          (*f)(t,y,n,d);
          for(i=0;i<=n-1;i++)
          {
              b=y;
              e=y;
          }
          for(k=0;k<=2;k++)
          {
              for(i=0;i<=n-1;i++)
              {
                  y=e+a[k]*d;
                  b=b+a[k+1]*d/3.0;
              }
              tt=t+a[k];
              (*f)(tt,y,n,d);
          }
          for(i=0; i<=n-1; i++)
              y=b+hh*d/6.0;
          t=t+dt;
      }
      p=0.0;
      for(i=0; i<=n-1;i++)
      {
          q=fabs(y-g);
          if(q>p) p=q;
      }
      hh=hh/2.0;
      m=m+m;
  }
  free(g);
  free(b);
  free(c);
  free(d);
  free(e);
  return;
}
**********************************************************************************************************
#include "stdio.h"
#include "math.h"
#include "stdlib.h"


#define PI 3.1415
#define g 9.8
int num_ball, n_rotor;//num_ball 滚珠的数量，n_rotor 转子转速

double damp,e,kn,w_rotor,w_c,F_out,m,gama,BN,T_vc,dt;
//damp 阻尼；e-偏心；kn-刚度；w_rotor 转子角速度；
//w_c 内圈转速；F_out-外力；m-转子质量；gama-游隙；
//BN,轴承参数，用以表示变刚度频率和转子转速的关系
//T_vc 变刚度的周期，60/(n_rotor*BN);%
//dt 周期内积分步长，一般一个周期取200个点
main()
{
  int i,j;
  void rkt2f(double,double [],int, double[]);
  double t,h,eps,y[4];

  damp=200.0;
  e=0.0;
  kn=7.055e9;
  n_rotor=10000;
  num_ball = 9;
  F_out=6.0;
  m=1.0;
  gama=2e-5;
  BN=3.6;
  w_rotor=n_rotor*PI/30.0;        //转速转换为角速度
  w_c=0.4*w_rotor;              //滚珠公转角速度
  T_vc=50.0/(3.0*n_rotor);        //变刚度周期T_vc=60/(n_rotor*3.6);
  dt=1.0/(12*n_rotor);             //一个周期取200个点0.005*T_vc=1/(12*n_rotor)
  
  t=0.0;
  h=dt;
  printf("dt=",T_vc);
  eps=1e-7;
  y[0]=0.0;y[1]=1e-8;y[2]=0.0;y[3]=1e-8;         //initial
  printf("\n");
  printf("t=%7.3f y[0]=%e y[1]=%e y[2]=%e y[3]=%e\n",t,y[0],y[1],y[2],y[3]);

  {
      for(i=1;i<=200000;i++)//计算1000个周期         
      { rkt2(t,h,y,4,eps,rkt2f);
       // mrsn(t,h,y,4,eps,rkt2f);
      t=t+h;      
      }
      printf("t=%7.3f y[0]=%e y[1]=%e y[2]=%e y[3]=%e\n", t,y[0],y[1],y[2],y[3]);
  }  
}

void rkt2f(t,y,n,d)
int n;
double t,y[],d[];
{
  int i,j,k;
  double *f;
  double *sita;                //滚珠的角度
  double *clearance;           //滚珠与内圈的间隙
  t=t;
  n=n;
  f=(double*)malloc(2*sizeof(double));
  sita=(double*)malloc(num_ball*sizeof(double));
  clearance=(double*)malloc(num_ball*sizeof(double));
  for(j=0;j<=1;j++)
      f[j]=0.0;
  for(k=0;k<num_ball;k++)
  {
      sita[k]=0.0;
      clearance[k]=0.0;
  }
  for(i=0;i<num_ball;i++)
  {
      sita=2*PI*(i-1)/num_ball+w_c*t;
      clearance=y[2]*sin(sita)+y[0]*cos(sita)-gama;
      if(clearance<=0)
          clearance=0;
      f[0]=f[0]+kn*pow(clearance,1.5)*cos(sita);
      f[1]=f[1]+kn*pow(clearance,1.5)*sin(sita);
  }
   
  d[0]=y[1];
  d[1]=e*w_rotor*w_rotor*cos(w_rotor*t)+g+F_out-damp*y[1]-f[0]-f[1];
  d[2]=y[3];
  d[3]=e*w_rotor*w_rotor*sin(w_rotor*t)-damp*y[3]-f[0]-f[1];
  free(f);
  free(sita);
  free(clearance);
  return;
}

胡晓宇

最后一个for循环中，sita为sita；clearance应为clearance，不知道怎么回事，显示的不正确，望谅解！

胡晓宇

最有一个for 循环中,sita 应为sita(sita数组的第i 个元素) clearance应为clearance(clearance数组的第i个元素)

sizhiyuan2006

谢谢分享