Я пытаюсь реализовать метод Эйлера, runge-kutta и метод средней точки в numpy Python.
Я хочу реализовать функцию integrateall, которая в зависимости от метода (euler, runge-kutta или middle-point) интегрируется по всем интервалам. Другие функции (euler, rk4 и middlepoint) должны сделать только один шаг и вернуть значения.
Вот мой код Python:
import numpy as np
import matplotlib as plt
from scipy.integrate import odeint
def oneStepProcess(y, t, h, ode):
'''
just to represents the form that one-step-processes that follow (euler, rk4 and middlepoint) should have
:param y: last calculated value at time t; y_k hast n components
y is an np.array
:param t: time t of last calculated step
:param h: timestep
:param ode: ode to solve as function
:return: y_new Solution at time t_{new}=t+h (as np.array)
'''
def euler(y, t, h, ode):
y_new = y + ode(t, y)*h
return y_new
def f1(t, y):
return -5*y
def fn(t, y):
dy = np.zeros_like(y)
dy[0] = -5*y[1]
dy[1] = 3*y[0]+y[1]
return dy
assert np.allclose(euler(np.array([1.]), 1., 0.1, f1), np.array([0.5]))
assert np.allclose(euler(np.array([1., -1]), 1., 0.2, fn), np.array([2, -0.6]))
def middlepoint(y, t, h, ode):
mid_step = 0.5*h
y_mid = y + (mid_step*ode(t, y))
ynew = y + h*ode(y_mid, t+mid_step)
return ynew
def rk4(y, t, h, ode):
mid_step = 0.5*h
k1 = y
k2 = ode(t+mid_step, y+(mid_step*k1))
k3 = ode(t+mid_step, y+(mid_step*k2))
k4 = ode(t+h, y+(h*k3))
ynew = y + (h*((k1/6)+(k2/3)+(k3/3)+(k4/6)))
return ynew
def integrateall(method, ode, y0, t0, tend, N, intermediate=False):
"""
:param method: one-step-process
:param ode: RHS of ode
:param y0: Start value
:param t0: Start time
:param tend: End time
:param N: Number of steps
:param intermediate: True, if we want to show intermediate solutions when plotting
:return: Solution at endtime (or solution AND intermediate times, if intermediate=True)
"""
tvals = [t0]
yvals = [y0]
t = t0
y = y0
h = (tend-t0/N)
while t < tend:
h = min(h, tend-t)
#t ,y = method(yvals[-1],tvals[-1],h,ode)
t ,y = method(y,t,h,ode)
tvals.append(t)
yvals.append(y)
yend = yvals[-1]
print(yvals)
if intermediate:
return np.array(yvals), np.array(tvals)
else:
return yend
def mouse(t, y):
dy = np.array([1 - (y[0]+y[1]), y[0]-y[1]])
return dy
def fsimple(t, y):
return t, y
y= integrateall(euler, fsimple, np.array([1.,]), 0., 1., 100, False)
# wanted to test the numpy-functions 'euler', 'rk4' and 'middlepoint' in 'integrateall' with the functions fsimple, mouse, f1 and fn but none works. Just left the first one (fsimple) for clarity. If you could help me why euler, rk4 and middlepoint are not working with 'integrateall' as well, it would be perfect!
Вот ошибка, которую я получаю:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-17-6624311b1e40> in <module>()
90 return t, y
91
---> 92 y= integrateall(euler, fsimple, np.array([1.,]), 0., 1., 100, False)
93
<ipython-input-17-6624311b1e40> in integrateall(method, ode, y0, t0, tend, N, intermediate)
71 h = min(h, tend-t)
72 #t ,y = method(yvals[-1],tvals[-1],h,ode)
---> 73 t ,y = method(y,t,h,ode)
74 tvals.append(t)
75 yvals.append(y)
<ipython-input-17-6624311b1e40> in euler(y, t, h, ode)
17
18 def euler(y, t, h, ode):
---> 19 y_new = y + ode(t, y)*h
20 return y_new
21
TypeError: can't multiply sequence by non-int of type 'float'