X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fcomputational_physics.git;a=blobdiff_plain;f=rk.c;fp=rk.c;h=65723287f180350c023e51c5a311e88b1c103cd8;hp=0000000000000000000000000000000000000000;hb=33422521a19e7afa578c599705bbfb0c8d31b88a;hpb=a649f7232da963470fd43e9a1a5384944ba5ed62

diff --git a/rk.c b/rk.c
new file mode 100644
index 0000000..6572328
--- /dev/null
+++ b/rk.c
@@ -0,0 +1,59 @@
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <math.h>
+
+#define N 10000
+#define M 100
+#define INT (2*M_PI)
+
+
+double k[4][2];
+double y[N+1][2];
+double tau;
+int i;
+
+double f0(double t,double y0,double y1) {
+  return y1;
+}
+
+double f1(double t,double y0,double y1) {
+  double q=.5,b=1.15,w=1.0*2/3;
+
+  return (-1.0*sin(y0)+b*cos(w*t)-q*y1);
+  // return -y0;
+}
+
+int main() {
+
+  y[0][0]=0;
+  y[0][1]=2;
+
+  tau=INT/M;
+
+  for(i=0;i<N;i++) {
+    printf("#run %d\n",i);
+
+    k[0][0]=tau*f0(i*tau,y[i][0],y[i][1]);
+    k[0][1]=tau*f1(i*tau,y[i][0],y[i][1]);
+
+    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]);
+    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]);
+
+    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]);
+    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]);
+
+    k[3][0]=tau*f0(i*tau+tau,y[i][0]+k[2][0],y[i][1]+k[2][1]);
+    k[3][1]=tau*f1(i*tau+tau,y[i][0]+k[2][0],y[i][1]+k[2][1]);
+
+    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);
+    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);
+
+  }
+
+  for(i=0;i<N;i++) printf("%f %f\n",y[i][0],y[i][1]);
+
+  return 1;
+}
+
+
+