rk.c

   1 #define _GNU_SOURCE
   2 #include <stdio.h>
   3 #include <math.h>
   4
   5 #define N 10000
   6 #define M 100
   7 #define INT (2*M_PI)
   8
   9
  10 double k[4][2];
  11 double y[N+1][2];
  12 double tau;
  13 int i;
  14
  15 double f0(double t,double y0,double y1) {
  16   return y1;
  17 }
  18
  19 double f1(double t,double y0,double y1) {
  20   double q=.5,b=1.15,w=1.0*2/3;
  21
  22   return (-1.0*sin(y0)+b*cos(w*t)-q*y1);
  23   // return -y0;
  24 }
  25
  26 int main() {
  27
  28   y[0][0]=0;
  29   y[0][1]=2;
  30
  31   tau=INT/M;
  32
  33   for(i=0;i<N;i++) {
  34     printf("#run %d\n",i);
  35
  36     k[0][0]=tau*f0(i*tau,y[i][0],y[i][1]);
  37     k[0][1]=tau*f1(i*tau,y[i][0],y[i][1]);
  38
  39     k[1][0]=tau*f0(i*tau+0.5*tau,y[i][0]+0.5*k[0][0],y[i][1]+0.5*k[0][1]);
  40     k[1][1]=tau*f1(i*tau+0.5*tau,y[i][0]+0.5*k[0][0],y[i][1]+0.5*k[0][1]);
  41
  42     k[2][0]=tau*f0(i*tau+0.5*tau,y[i][0]+0.5*k[1][0],y[i][1]+0.5*k[1][1]);
  43     k[2][1]=tau*f1(i*tau+0.5*tau,y[i][0]+0.5*k[1][0],y[i][1]+0.5*k[1][1]);
  44
  45     k[3][0]=tau*f0(i*tau+tau,y[i][0]+k[2][0],y[i][1]+k[2][1]);
  46     k[3][1]=tau*f1(i*tau+tau,y[i][0]+k[2][0],y[i][1]+k[2][1]);
  47
  48     y[i+1][0]=(y[i][0]+(k[0][0]+2*k[1][0]+2*k[2][0]+k[3][0])/6); //%(2*M_PI);
  49     y[i+1][1]=(y[i][1]+(k[0][1]+2*k[1][1]+2*k[2][1]+k[3][1])/6); //%(2*M_PI);
  50
  51   }
  52
  53   for(i=0;i<N;i++) printf("%f %f\n",y[i][0],y[i][1]);
  54
  55   return 1;
  56 }
  57
  58
  59