one-file-projects/heun.py

30 lines
650 B
Python
Raw Normal View History

2015-05-23 23:44:26 +02:00
# Verfahren von Heun zur numerischen Lösung von Anfangswertaufgaben.
def run_heun(K,h,x0,y0,ys):
if not isinstance(y0,list): y0 = [y0]
2015-05-26 19:46:44 +02:00
if not isinstance(ys,list): ys = [ys]
2015-05-23 23:44:26 +02:00
xk = x0
yk = y0
2015-05-26 19:46:44 +02:00
xka = 0
yka = [0]*len(yk)
ykp = [0]*len(yk)
xks = [x0]
2015-05-23 23:44:26 +02:00
yks = [ [v] for v in y0 ]
for i in range(K):
2015-05-26 19:46:44 +02:00
xka = xk + h
xks.append(xka)
for j in range(len(ys)):
yka[j] = yk[j] + h * ys[j](xk,*yk)
2015-05-23 23:44:26 +02:00
for j in range(len(ys)):
2015-05-26 19:46:44 +02:00
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
2015-05-23 23:44:26 +02:00