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 +