From 0396c93cad600949ae93ef2ce6d6e024c3bd0778 Mon Sep 17 00:00:00 2001
From: Valentin Gehrke <madmaurice@zom.bi>
Date: Tue, 26 May 2015 19:46:44 +0200
Subject: [PATCH] =?UTF-8?q?Reimplementierung=20f=C3=BCr=20DGlS?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 heun.py | 50 +++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 39 insertions(+), 11 deletions(-)

diff --git a/heun.py b/heun.py
index 262dd60..f36bfe9 100644
--- a/heun.py
+++ b/heun.py
@@ -1,22 +1,50 @@
 # Verfahren von Heun zur numerischen Lösung von Anfangswertaufgaben.
-
-def heun(xk,yk,h,ys):
-    xk1 = xk + h
-    yk1p = yk + h * ys(xk,yk)
-    yk1 = 0.5*(yk+yk1p+h*ys(xk1,yk1p))
-    return (xk1,yk1p)
+#
+#def heun(xk,yk,h,ys):
+#    xk1 = xk + h
+#    yk1p = yk + h * ys(xk,yk)
+#    yk1 = 0.5*(yk+yk1p+h*ys(xk1,yk1p))
+#    return (xk1,yk1p)
+#
+#def run_heun(K,h,x0,y0,ys):
+#    if not isinstance(ys,list): ys = [ys]
+#    if not isinstance(y0,list): y0 = [y0]
+#    xk = x0
+#    yk = y0
+#    xks = x0
+#    yks = [ [v] for v in y0 ]
+#
+#    for i in range(K):
+#        for j in range(len(ys)):
+#            xk,yk[j] = heun(xk,yk[j],h,ys[i])
+#            xks.append(xk)
+#            yks[j].append(yk[j])
+#
 
 def run_heun(K,h,x0,y0,ys):
-    if not isinstance(ys,list): ys = [ys]
     if not isinstance(y0,list): y0 = [y0]
+    if not isinstance(ys,list): ys = [ys]
+
     xk = x0
     yk = y0
-    xks = x0
+    xka = 0
+    yka = [0]*len(yk)
+    ykp = [0]*len(yk)
+
+    xks = [x0]
     yks = [ [v] for v in y0 ]
 
     for i in range(K):
+        xka = xk + h
+        xks.append(xka)
         for j in range(len(ys)):
-            xk,yk[j] = heun(xk,yk[j],h,ys[i])
-            xks.append(xk)
-            yks[j].append(yk[j])
+            yka[j] = yk[j] + h * ys[j](xk,*yk)
+        for j in range(len(ys)):
+            ykp[j] = 0.5*(yk[j] + yka[j] + h * ys[j](xka,*yka))
+            yks[j].append(ykp[j])
+        yk = ykp
+        xk = xka
+
+    return [xks] + yks
+