|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 |
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]) |
跳至内容
Personal Blog
浙ICP备2021019730-1 
浙公网安备 33010902002953号
Copyright © 2024 PanCake