added runge kutta example

diff --git a/Makefile b/Makefile
index 0af42366cdb042edf95631ff7ee6d5d6440de76c..e08fef7382438945c5443f885777cc19f9b18500 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,11 +1,12 @@
 # computational physics Makefile, rules for all example programs
 
+CC=gcc
 INCLUDEDIR = /usr/include
 CFLAGS = -O3 -Wall
 LIBS = -L/usr/lib -lm
 
 API = g_plot.o general.o
-OBJS = newton zentral homogen integral-1_2 integral-2_2 polynom_interpolation kettenbruchentwicklung bessel_1 bessel_2 check_rand mc_int
+OBJS = newton zentral homogen integral-1_2 integral-2_2 polynom_interpolation kettenbruchentwicklung bessel_1 bessel_2 check_rand mc_int nullstellen rk
 
 all: $(OBJS)
 
@@ -42,7 +43,13 @@ check_rand: $(API)
 mc_int: $(API)
        $(CC) $(CFLAGS) -o $@ $(API) $(LIBS) mc_int.c
 
+nullstellen: $(API)
+       $(CC) $(CFLAGS) -o $@ $(API) $(LIBS) nullstellen.c
+
+rk: $(API)
+       $(CC) $(CFLAGS) -o $@ $(API) $(LIBS) rk.c
+
 clean:
-       rm $(API) $(OBJS)
+       rm -f $(API) $(OBJS)
 
 remake: clean all

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;
+}
+
+
+