全部是正确答案，慎用。Adams和Milne需要用另外的算法如runge-kutta（代码示例求法）求前3个数，当然也可以用Euler和Euler Improved，大家不要都一样啦。Python缩进规范请注意！
Python Euler Improved
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
def solveByEulerImproved(f, epsilon, a, y_a, b): function = Result.get_function(f) n = 1 delta = float("inf") while abs(delta) >= epsilon / 10: n *= 2 h = (b - a) / n pre_y = y_a pre_fy = function(a, y_a) for i in range(n): x = a + (i + 1 / 2) * h yHalf = pre_y + pre_fy * h / 2 fyHalf = function(x, yHalf) y = pre_y + h * fyHalf delta = pre_y - y pre_y = y pre_fy = function(x + h / 2, y) return y |
Python Euler
Euler method 理论上代码正确，但是实际上因为算法的原因，与测试样例的值偏差会过大，要通过测试样例只需要初始化时把 n = 10000000
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
def solveByEuler(f, epsilon, a, y_a, b): function = Result.get_function(f) n = 1 delta = float("inf") while abs(delta) >= epsilon / 10: n *= 2 h = (b - a) / n pre_y = y_a for i in range(n): x = a + i * h fy = function(x, pre_y) y = pre_y + h * fy delta = pre_y - y pre_y = y return y |
Python Runge-Kutta
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
def solveByRungeKutta(f, epsilon, a, y_a, b): function = Result.get_function(f) n = 1 delta = float("inf") while abs(delta) >= epsilon / 10: n *= 2 h = (b - a) / n pre_y = y_a for i in range(n): x = a + i * h k1 = h * function(x, pre_y) k2 = h * function(x + h / 2, pre_y + k1 / 2) k3 = h * function(x + h / 2, pre_y + k2 / 2) k4 = h * function(x + h, pre_y + k3) y = (k1 + 2 * k2 + 2 * k3 + k4) / 6 + pre_y delta = pre_y - y pre_y = y return y |
Python Adams
|
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 |
def runge_kutta(h, a, y_a, function): y_dot = [function(a, y_a)] pre_y = y_a for i in range(3): x = a + i * h k1 = h * function(x, pre_y) k2 = h * function(x + h / 2, pre_y + k1 / 2) k3 = h * function(x + h / 2, pre_y + k2 / 2) k4 = h * function(x + h, pre_y + k3) y = (k1 + 2 * k2 + 2 * k3 + k4) / 6 + pre_y y_dot.append(function(x + h, y)) pre_y = y return y, y_dot def solveByAdams(f, epsilon, a, y_a, b): function = Result.get_function(f) n = 2 delta = float("inf") while abs(delta) >= epsilon / 10: n *= 2 h = (b - a) / n y, y_dot = Result.runge_kutta(h, a, y_a, function) for i in range(n - 3): delta = (55 * y_dot[-1] - 59 * y_dot[-2] + 37 * y_dot[-3] - 9 * y_dot[-4]) / 24 * h y += delta y_dot.append(function(a + (i + 4) * h, y)) return y |
Python Milne
|
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 |
def runge_kutta(h, a, y_a, function): y_dot = [function(a, y_a)] y_ = [y_a] pre_y = y_a for i in range(3): x = a + i * h k1 = h * function(x, pre_y) k2 = h * function(x + h / 2, pre_y + k1 / 2) k3 = h * function(x + h / 2, pre_y + k2 / 2) k4 = h * function(x + h, pre_y + k3) y = (k1 + 2 * k2 + 2 * k3 + k4) / 6 + pre_y y_dot.append(function(x + h, y)) y_.append(y) pre_y = y return y_, y_dot def solveByMilne(f, epsilon, a, y_a, b): function = Result.get_function(f) n = 2 delta = float("inf") while abs(delta) >= epsilon / 10: n *= 2 h = (b - a) / n y, y_dot = Result.runge_kutta(h, a, y_a, function) for i in range(n - 3): y.append(y[-4] + 4 * h / 3 * (2 * y_dot[-1] - y_dot[-2] + 2 * y_dot[-3])) y_dot.append(function(a + (i + 4) * h, y[-1])) y[-1] = y[-3] + h / 3 * (y_dot[-1] + 4 * y_dot[-2] + y_dot[-3]) delta = y[-1] - y[-2] y_dot[-1] = function(a + (i + 4) * h, y[-1]) return y[-1] |
浙公网安备 33010902002953号