求积分以及方程的解需要自己修改f(x)函数里面的返回值。求非线性方程的解需要安装sympy包。
牛顿法解一元非线性方程
|
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 |
from math import * from sympy import * def f(x): return x ** 2 - 3 * x ** 3 def newtonByStep(x0, step): x = symbols("x", real=True) _f = f(x) f_x = diff(_f, x) for i in range(step): x1 = x0 - f(x0) / f_x.subs({x: x0}) print("Step{}: x = {}, f(x) = {}".format(i + 1, x1, f(x1))) x0 = x1 def newtonByEpsilon(x0, epsilon): x = symbols("x", real=True) _f = f(x) f_x = diff(_f, x) error = float("inf") while error >= epsilon: x1 = x0 - f(x0) / f_x.subs({x: x0}) error = abs(x1 - x0) x0 = x1 return x0 if __name__ == "__main__": x0 = 1.04 newtonByStep(x0, 5) print(newtonByEpsilon(x0, 0.001)) |
牛顿法解二元非线性方程
|
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 |
from math import * from sympy import * def f(x, y): return 2 * x ** 3 - y ** 2 - 1 def g(x, y): return x * y ** 3 - y - 4 def newtonByStep(x0, y0, step): x, y = symbols("x y", real=True) _f = f(x, y) _g = g(x, y) f_x = diff(_f, x) f_y = diff(_f, y) g_x = diff(_g, x) g_y = diff(_g, y) for i in range(step): x1 = x0 + (f(x0, y0) * g_y.subs({x: x0, y: y0}) - g(x0, y0) * f_y.subs({x: x0, y: y0})) / ( g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0})) y1 = y0 + (g(x0, y0) * f_x.subs({x: x0, y: y0}) - f(x0, y0) * g_x.subs({x: x0, y: y0})) / ( g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0})) print("Step{}: x = {}, y = {}, f(x, y) = {}, g(x, y) = {}".format(i + 1, x1, y1, f(x1, y1), g(x1, y1))) x0 = x1 y0 = y1 def newtonByEpsilon(x0, y0, epsilon): x, y = symbols("x y", real=True) _f = f(x, y) _g = g(x, y) f_x = diff(_f, x) f_y = diff(_f, y) g_x = diff(_g, x) g_y = diff(_g, y) error = float("inf") while error >= epsilon: x1 = x0 + (f(x0, y0) * g_y.subs({x: x0, y: y0}) - g(x0, y0) * f_y.subs({x: x0, y: y0})) / ( g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0})) y1 = y0 + (g(x0, y0) * f_x.subs({x: x0, y: y0}) - f(x0, y0) * g_x.subs({x: x0, y: y0})) / ( g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0})) error = max(abs(x1 - x0), abs(y1 - y0)) x0 = x1 y0 = y1 return x0, y0 if __name__ == "__main__": x0 = 1.2 y0 = 1.7 newtonByStep(x0, y0, 5) print(newtonByEpsilon(x0, y0, 0.001)) |
5种积分方法
|
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 |
from math import * def sgn(x): if x > 0: return 1 elif x == 0: return 0 else: return -1 def func(x): # return sgn(sin(pi / sin(x))) return x ** 2 + 2 def simpson(begin, between, i): a = begin + (i - 1) * between b = begin + i * between return between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6 def rightRectangle(begin, between, i): return between * func(begin + between * i) def leftRectangle(begin, between, i): return between * func(begin + between * (i - 1)) def midRectangle(begin, between, i): return between * func(begin + between / 2 * (2 * i - 1)) def trapezoidRectangle(begin, between, i): return between / 2 * (func(begin + between * (i - 1)) + func(begin + between * i)) def cal_IntegralByEpsilon(begin, end, epsilon, method): n = 1 result = 0 preResult = float("inf") while abs(preResult - result) >= epsilon: preResult = result result = 0 n *= 2 between = (end - begin) / n for i in range(n): try: result += method(begin, between, i + 1) except: return "Integrated function has discontinuity or does not defined in current interval" return result def cal_IntegralByN(begin, end, n, method): result = 0 between = (end - begin) / n for i in range(n): try: result += method(begin, between, i + 1) except: return "Integrated function has discontinuity or does not defined in current interval" if __name__ == "__main__": begin = 0.2 end = 0.5 epsilon = 0.0001 print(cal_IntegralByEpsilon(begin, end, epsilon, simpson)) print(cal_IntegralByEpsilon(begin, end, epsilon, rightRectangle)) print(cal_IntegralByEpsilon(begin, end, epsilon, leftRectangle)) print(cal_IntegralByEpsilon(begin, end, epsilon, midRectangle)) print(cal_IntegralByEpsilon(begin, end, epsilon, trapezoidRectangle)) begin = 0.2 end = 0.5 n = 10 print(cal_IntegralByN(begin, end, n, simpson)) print(cal_IntegralByN(begin, end, n, rightRectangle)) print(cal_IntegralByN(begin, end, n, leftRectangle)) print(cal_IntegralByN(begin, end, n, midRectangle)) print(cal_IntegralByN(begin, end, n, trapezoidRectangle)) |
绝对误差与相对误差
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
def errorByreal(real, standard): absolute = abs(real - standard) realative = absolute / real print("Absolute error = %e \nRelative error = %e" % (absolute, realative)) def errorByEpsilon(real, epsilon): accuracy = len(str(epsilon).split(".")[-1]) standard = round(real, accuracy) absolute = abs(real - standard) realative = absolute / real print("Absolute error = %e \nRelative error = %e" % (absolute, realative)) if __name__ == "__main__": real = 0.63840532 epsilon = 0.01 errorByEpsilon(real, epsilon) real = 0.63840532 standard = 0.6384 errorByreal(real, standard) |
最小二乘法一元多项式拟合
|
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 |
def gauss(matrix): array = [] for i in range(0, len(matrix) - 1): aii = matrix[i][i] if aii == 0: continue for j in range(i + 1, len(matrix)): aji = matrix[j][i] for k in range(i, len(matrix[j])): matrix[j][k] = matrix[j][k] - aji / aii * matrix[i][k] for i in range(len(matrix) - 1, -1, -1): if matrix[i][i] == 0: exit("The system has no roots of equations or has an infinite set of them.") array.append(matrix[i][-1] / matrix[i][i]) for j in range(0, i): matrix[j][-1] -= matrix[j][i] * array[-1] array.reverse() return array def leastSquareRegression(points, order): result = [] for i in range(order + 1): sum = 0 for j in range(len(points)): sum += points[j][0] ** i * points[j][1] result.append(sum) matrix = [] for i in range(order + 1): matrix.append([]) for j in range(2 * order): sum = 0 for k in range(len(points)): sum += points[k][0] ** (i + j) matrix[i].append(sum) matrix[i].append(result[i]) return gauss(matrix) def R2_func(y_test, y): return 1 - ((y_test - y) ** 2).sum() / ((y.mean() - y) ** 2).sum() def test(points, result): import numpy as np approach = [] for i in range(len(points)): sum = 0 for j in range(order + 1): sum += points[i][0] ** j * result[j] approach.append(sum) print("R^2 result =", R2_func(np.array(approach), np.array([i[1] for i in points]))) if __name__ == "__main__": points = [[1, 2], [3, 4], [5, 77], [8, -3]] order = 2 result = leastSquareRegression(points, order) print("y =", " + ".join(["{}X^{}".format(result[i], i) for i in range(order + 1)])) test(points, result) |
浙公网安备 33010902002953号