added runge kutta example

[physik/computational_physics.git] / rk.c

diff --git a/rk.c b/rk.c
new file mode 100644 (file)
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;
+}
+
+
+