class ODE:
def __init__(self, function, x0, y0):
self.function = function
self.x0 = x0
self.y0 = y0
def euler(self, h, step):
result = []
pre_y = self.y0
for i in range(step):
x = self.x0 + i * h
fy = self.function(x, pre_y)
y = pre_y + h * fy
result.append(y)
pre_y = y
# print("y_{} = {}".format(i + 1, y))
return result
def eulerImproved(self, h, step):
result = []
pre_y = self.y0
pre_fy = self.y0
for i in range(step):
x = self.x0 + (i + 1 / 2) * h
yHalf = pre_y + pre_fy * h / 2
fyHalf = self.function(x, yHalf)
y = pre_y + h * fyHalf
result.append(y)
pre_y = y
pre_fy = function(x + h / 2, y)
# print("y_{} = {}".format(i + 1, y))
return result
def runge_kutta(self, h, step):
result = []
pre_y = self.y0
for i in range(step):
x = self.x0 + i * h
k1 = h * self.function(x, pre_y)
k2 = h * self.function(x + h / 2, pre_y + k1 / 2)
k3 = h * self.function(x + h / 2, pre_y + k2 / 2)
k4 = h * self.function(x + h, pre_y + k3)
y = (k1 + 2 * k2 + 2 * k3 + k4) / 6 + pre_y
result.append(y)
pre_y = y
# print("y_{} = {}".format(i + 1, y))
return result
def milne(self, h, step):
y = [self.y0] + self.runge_kutta(h, 3)
y_dot = [function(self.x0 + i * h, y[i]) for i in range(4)]
for i in range(step - 3):
y.append(y[-4] + 4 * h / 3 * (2 * y_dot[-1] - y_dot[-2] + 2 * y_dot[-3]))
y_dot.append(function(self.x0 + (i + 4) * h, y[-1]))
y[-1] = y[-3] + h / 3 * (y_dot[-1] + 4 * y_dot[-2] + y_dot[-3])
y_dot[-1] = function(self.x0 + (i + 4) * h, y[-1])
# print("y_{} = {}".format(i + 4, y[-1]))
return y[1:]
def adams(self, h, step):
y = [self.y0] + self.runge_kutta(h, 3)
y_dot = [function(self.x0 + i * h, y[i]) for i in range(4)]
for i in range(step - 3):
y.append(y[-1] + (55 * y_dot[-1] - 59 * y_dot[-2] + 37 * y_dot[-3] - 9 * y_dot[-4]) / 24 * h)
y_dot.append(function(self.x0 + (i + 4) * h, y[-1]))
# print("y_{} = {}".format(i + 4, y[-1]))
return y[1:]
if __name__ == "__main__":
function = lambda x, y: y - 2 * x / y
x0 = 0
y0 = 1
h = 0.002
step = 1000
ode = ODE(function, x0, y0)
print(ode.euler(h, step)[-1])
print("--------------------------------------------")
print(ode.eulerImproved(h, step)[-1])
print("--------------------------------------------")
print(ode.runge_kutta(h, step)[-1])
print("--------------------------------------------")
print(ode.adams(h, step)[-1])
print("--------------------------------------------")
print(ode.milne(h, step)[-1])