24 Hours – 第 7 页 – Pancake's Personal Website

Python 牛顿法求解多元非线性方程组

代码使用方式可见代码中示例代码。需安装numpy和sympy库。

from math import *
from sympy import *
import sympy
import numpy as np


class Newton:
    def __init__(self, equations):
        self.equations = equations
        self.n = len(equations)
        self.x = sympy.symbols(" ".join("x{}".format(i) for i in range(self.n)) + " x{}".format(self.n), real=True)
        self.equationsSymbol = [equations[i](self.x) for i in range(self.n)]
        # 初始化 Jacobian 矩阵
        self.J = np.zeros(self.n * self.n, dtype=sympy.core.add.Add).reshape(self.n, self.n)
        for i in range(self.n):
            for j in range(self.n):
                self.J[i][j] = sympy.diff(self.equationsSymbol[i], self.x[j])

    def cal_J(self, x):
        dict = {self.x[i]: x[i] for i in range(self.n)}
        J = np.zeros(self.n * self.n).reshape(self.n, self.n)
        for i in range(self.n):
            for j in range(self.n):
                J[i][j] = self.J[i][j].subs(dict)
        return J

    def cal_f(self, x):
        f = np.zeros(self.n)
        for i in range(self.n):
            f[i] = self.equations[i](x)
        f.reshape(self.n, 1)
        return f

    def newtonByStep(self, x0, step):
        x0 = np.array(x0)
        for i in range(step):
            x0 = x0 - np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)
            print("Step {}:".format(i + 1), ", ".join(["x{} = {}".format(j + 1, x0[j]) for j in range(self.n)]))
        return x0

    def newtonByEpsilon(self, x0, epsilon):
        error = float("inf")
        while error >= epsilon:
            cal = np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)
            error = max(abs(cal))
            x0 = x0 - cal
        print(x0)
        return x0


if __name__ == "__main__":
    # equations 为 equation = 0 的方程组
    # x0 为初始值 step 为迭代次数 epsilon 为精度

    # 多元非线性使用方法
    equations = [lambda x: cos(0.4 * x[1] + x[0] ** 2) + x[0] ** 2 + x[1] ** 2 - 1.6,
                 lambda x: 1.5 * x[0] ** 2 - x[1] ** 2 / 0.36 - 1,
                 lambda x: 3 * x[0] + 4 * x[1] + 5 * x[2]]

    newton = Newton(equations)
    newton.newtonByStep([1, 1, 1], 5)
    newton.newtonByEpsilon([1, 1, 1], 0.001)

    # 一元非线性使用方法
    equations = [lambda x: cos(x[0]) + sin(x[0])]

    newton = Newton(equations)
    newton.newtonByStep([1], 5)
    newton.newtonByEpsilon([1], 0.001)

    # 考试题
    equations = [lambda x: 2 * x[0] ** 3 - x[1] ** 2 - 1,
                 lambda x: x[0] * x[1] ** 3 - x[1] - 4]

    newton = Newton(equations)
    newton.newtonByStep([1.2, 1.7], 5)

from math import *

from sympy import *

import sympy

import numpy as np

class Newton:

def __init__(self, equations):

self.equations = equations

self.n = len(equations)

self.x = sympy.symbols(" ".join("x{}".format(i) for i in range(self.n)) + " x{}".format(self.n), real=True)

self.equationsSymbol = [equations[i](self.x) for i in range(self.n)]

# 初始化 Jacobian 矩阵

self.J = np.zeros(self.n * self.n, dtype=sympy.core.add.Add).reshape(self.n, self.n)

for i in range(self.n):

for j in range(self.n):

self.J[i][j] = sympy.diff(self.equationsSymbol[i], self.x[j])

def cal_J(self, x):

dict = {self.x[i]: x[i] for i in range(self.n)}

J = np.zeros(self.n * self.n).reshape(self.n, self.n)

for i in range(self.n):

for j in range(self.n):

J[i][j] = self.J[i][j].subs(dict)

return J

def cal_f(self, x):

f = np.zeros(self.n)

for i in range(self.n):

f[i] = self.equations[i](x)

f.reshape(self.n, 1)

return f

def newtonByStep(self, x0, step):

x0 = np.array(x0)

for i in range(step):

x0 = x0 - np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)

print("Step {}:".format(i + 1), ", ".join(["x{} = {}".format(j + 1, x0[j]) for j in range(self.n)]))

return x0

def newtonByEpsilon(self, x0, epsilon):

error = float("inf")

while error >= epsilon:

cal = np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)

error = max(abs(cal))

x0 = x0 - cal

print(x0)

return x0

if __name__ == "__main__":

# equations 为 equation = 0 的方程组

# x0 为初始值 step 为迭代次数 epsilon 为精度

# 多元非线性使用方法

equations = [lambda x: cos(0.4 * x[1] + x[0] ** 2) + x[0] ** 2 + x[1] ** 2 - 1.6,

lambda x: 1.5 * x[0] ** 2 - x[1] ** 2 / 0.36 - 1,

lambda x: 3 * x[0] + 4 * x[1] + 5 * x[2]]

newton = Newton(equations)

newton.newtonByStep([1, 1, 1], 5)

newton.newtonByEpsilon([1, 1, 1], 0.001)

# 一元非线性使用方法

equations = [lambda x: cos(x[0]) + sin(x[0])]

newton = Newton(equations)

newton.newtonByStep([1], 5)

newton.newtonByEpsilon([1], 0.001)

# 考试题

equations = [lambda x: 2 * x[0] ** 3 - x[1] ** 2 - 1,

lambda x: x[0] * x[1] ** 3 - x[1] - 4]

newton = Newton(equations)

newton.newtonByStep([1.2, 1.7], 5)

计算数学课程参考 Python代码

求积分以及方程的解需要自己修改f(x)函数里面的返回值。求非线性方程的解需要安装sympy包。

牛顿法解一元非线性方程

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))

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))

牛顿法解二元非线性方程

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))

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种积分方法

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))

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))

绝对误差与相对误差

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)

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)

最小二乘法一元多项式拟合

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)

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)

Python 求积分的几种方法

不同的方法最大区别在于公式的不同/矩形区域的不同，依据公式而定，需要什么代码你可以自己进行更改，很简单。
本文代码为不严格代码，最严格的代码为先判断连续性。当你需要对大量函数积分（懒得写函数连续性检测的时候）可以使用本文代码用try…except…不严格替代。
看不懂公式的可以去 CSDN 查看图片。

Simpson积分法:
\(\int _ { a } ^ { b } f ( x ) d x \approx \frac { ( b – a ) ( f ( a ) + f ( b ) + 4 f ( \frac { a + b } { 2 } ) ) } { 6 }\)

 def calculate_integral(begin, end, func, epsilon):
        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:
                    a = begin + i * between
                    b = begin + (i + 1) * between
                    result += between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6
                except:
                    return "Integrated function has discontinuity or does not defined in current interval"
        return result

def calculate_integral(begin, end, func, epsilon):

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:

a = begin + i * between

b = begin + (i + 1) * between

result += between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6

except:

return "Integrated function has discontinuity or does not defined in current interval"

return result

右矩形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { n } f ( a + \frac { b – a } { n } i )\)

def calculate_integral(begin, end, func, epsilon):
    n = 1
    result = 0
    preResult = float("inf")
    while abs(preResult - result) > epsilon:
        preResult = result
        result = 0
        n *= 2
        for i in range(n):
            try:
                result += (end - begin) / n * func(begin + (end - begin) / n * (i + 1))
            except:
                return "Integrated function has discontinuity or does not defined in current interval"
    return result

def calculate_integral(begin, end, func, epsilon):

n = 1

result = 0

preResult = float("inf")

while abs(preResult - result) > epsilon:

preResult = result

result = 0

n *= 2

for i in range(n):

try:

result += (end - begin) / n * func(begin + (end - begin) / n * (i + 1))

except:

return "Integrated function has discontinuity or does not defined in current interval"

return result

左矩形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { n } f ( a + \frac { b – a } { n } (i – 1) )\)

中矩形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { n } f ( a + \frac { (b – a)} { 2 n }(2i – 1 ) )\)

梯形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { 2n } (f ( a + \frac { b – a } { n }( i – 1) ) + f ( a + \frac { b – a } { n } i ))\)

Timus #1628. White Streaks

题目：有m*n的矩形，每个方块默认为白色，现给k个方块图上黑色，问有多少连续的num*1或1*num的长方形。如果一个1*1的方块周围没有白色方块，那么这个方块也包含在答案里面。

解析：行列遍历。对行和列分别进行查找，对行和列查找过程中单独成一个方块的坐标分别记录在两个集合中，最终只需要在结果上加上两个集合的交集中元素的个数。

#include<iostream>
#include<algorithm>
#include<iterator>
#include<set>
using namespace std;

struct Point
{
	int x, y;
	bool operator< (const Point & a)const
	{
		if (this->x == a.x)return this->y < a.y;
		return this->x < a.x;
	}
}point[60000];

int fx(Point& a, Point& b) { 
	if (a.x == b.x)return a.y < b.y;
	return a.x < b.x; 
}

int fy(Point& a, Point& b) {
	if (a.y == b.y)return a.x < b.x;
	return a.y < b.y;
}

void insert(set<Point>& a, int x, int y) {
	Point tempPoint;
	tempPoint.x = x;
	tempPoint.y = y;
	a.insert(tempPoint);

}

int main() {
	int m, n, k, ans = 0;
	cin >> n >> m >> k;
	set<Point> A;
	set<Point> B;
	set<Point> C;
	for (int i = 0; i < k; i++)cin >> point[i].x >> point[i].y;
	int now = 0, row = 0, col = 0;

	if (n == 1) {
		sort(point, point + k, fx);
		if (point[0].y == 1) ans++; 
		for (int i = 1; i < k; i++) {
			if (point[i].y - point[i - 1].y == 1) ans++;
		}
		if (point[k - 1].y == m) ans++;
		cout << k + 1 - ans;
		return 0;
	}

	if (m == 1) {
		sort(point, point + k, fy);
		if (point[0].x == 1) ans++;
		for (int i = 1; i < k; i++) {
			if (point[i].x - point[i - 1].x == 1) ans++;
		}
		if (point[k - 1].x == n) ans++;
		cout << k + 1 - ans;
		return 0;
	}

	sort(point, point + k, fx);
	while (row < n) {
		if (point[now].x != row + 1) {
			ans++;
			row++;
			continue;
		}
		if (point[now].y == 2) insert(A, point[now].x, 1);
		else if (point[now].y > 2) ans++;
		now++;
		while (now < k)
		{
			if (point[now].x == point[now - 1].x) {
				if (point[now].y - point[now - 1].y == 2) insert(A, point[now].x, point[now].y - 1);
				else if (point[now].y - point[now - 1].y > 2) ans++;
			}
			else {
				break;
			}
			now++;
		}
		if (point[now - 1].y == m - 1) insert(A, point[now - 1].x, m);
		else if (point[now - 1].y < m - 1) ans++;
		row++;
	}

	sort(point, point + k, fy);
	now = 0;
	while (col < m) {
		if (point[now].y != col + 1) {
			ans++;
			col++;
			continue;
		}
		if (point[now].x == 2) insert(B, 1, point[now].y);
		else if (point[now].x > 2) ans++;
		now++;
		while (now < k)
		{
			if (point[now].y == point[now - 1].y) {
				if (point[now].x - point[now - 1].x == 2) insert(B, point[now].x - 1, point[now].y);
				else if (point[now].x - point[now - 1].x > 2) ans++;
			}
			else {
				break;
			}
			now++;
		}
		if (point[now - 1].x == n - 1) insert(B, n, point[now - 1].y);
		else if (point[now - 1].x < n - 1) ans++;
		col++;
	}

	set_intersection(A.begin(), A.end(), B.begin(), B.end(), inserter(C, C.begin()));
	ans += C.size();
	cout << ans << endl;

}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

#include<iostream>

#include<algorithm>

#include<iterator>

#include<set>

using namespace std;

struct Point

{

int x, y;

bool operator< (const Point & a)const

{

if (this->x == a.x)return this->y < a.y;

return this->x < a.x;

}

}point[60000];

int fx(Point& a, Point& b) {

if (a.x == b.x)return a.y < b.y;

return a.x < b.x;

}

int fy(Point& a, Point& b) {

if (a.y == b.y)return a.x < b.x;

return a.y < b.y;

}

void insert(set<Point>& a, int x, int y) {

Point tempPoint;

tempPoint.x = x;

tempPoint.y = y;

a.insert(tempPoint);

}

int main() {

int m, n, k, ans = 0;

cin >> n >> m >> k;

set<Point> A;

set<Point> B;

set<Point> C;

for (int i = 0; i < k; i++)cin >> point[i].x >> point[i].y;

int now = 0, row = 0, col = 0;

if (n == 1) {

sort(point, point + k, fx);

if (point[0].y == 1) ans++;

for (int i = 1; i < k; i++) {

if (point[i].y - point[i - 1].y == 1) ans++;

}

if (point[k - 1].y == m) ans++;

cout << k + 1 - ans;

return 0;

}

if (m == 1) {

sort(point, point + k, fy);

if (point[0].x == 1) ans++;

for (int i = 1; i < k; i++) {

if (point[i].x - point[i - 1].x == 1) ans++;

}

if (point[k - 1].x == n) ans++;

cout << k + 1 - ans;

return 0;

}

sort(point, point + k, fx);

while (row < n) {

if (point[now].x != row + 1) {

ans++;

row++;

continue;

}

if (point[now].y == 2) insert(A, point[now].x, 1);

else if (point[now].y > 2) ans++;

now++;

while (now < k)

{

if (point[now].x == point[now - 1].x) {

if (point[now].y - point[now - 1].y == 2) insert(A, point[now].x, point[now].y - 1);

else if (point[now].y - point[now - 1].y > 2) ans++;

}

else {

break;

}

now++;

}

if (point[now - 1].y == m - 1) insert(A, point[now - 1].x, m);

else if (point[now - 1].y < m - 1) ans++;

row++;

}

sort(point, point + k, fy);

now = 0;

while (col < m) {

if (point[now].y != col + 1) {

ans++;

col++;

continue;

}

if (point[now].x == 2) insert(B, 1, point[now].y);

else if (point[now].x > 2) ans++;

now++;

while (now < k)

{

if (point[now].y == point[now - 1].y) {

if (point[now].x - point[now - 1].x == 2) insert(B, point[now].x - 1, point[now].y);

else if (point[now].x - point[now - 1].x > 2) ans++;

}

else {

break;

}

now++;

}

if (point[now - 1].x == n - 1) insert(B, n, point[now - 1].y);

else if (point[now - 1].x < n - 1) ans++;

col++;

}

set_intersection(A.begin(), A.end(), B.begin(), B.end(), inserter(C, C.begin()));

ans += C.size();

cout << ans << endl;

}

Timus #1067. Disk Tree

题目：输入N个文件/文件夹位置关系，输出层级关系。看题目例子就可以看懂。

解析：利用树状结构

代码缺点：该代码内存占用太多。在insert函数中，可以选择更高效的查找及插入算法，但这道题对运行时间要求没那没高，可以写的相对简单一点。Tree* child[70]这个数字70根据题目可以调整，只是越大消耗内存越大，可以存放的子节点越多，Timus这道题写70就够了。

代码优化建议：
1. 把Tree中改成vector<pair<string, Tree>> child;这样可以大大减少内存以及提供更高效的对string排序操作。
2. 把Tree中改成set<Tree> child;并对<运算符(根据string name;)进行重载。set提供了高效的insert和find函数，这种方法更高效。

#include<iostream>
#include<algorithm>
#include<vector>
using namespace std;

class Tree {
public:
	Tree() {
		level = 0;
	}
	Tree(const Tree& tree) {
		this->level = tree.level + 1;
	}
	int level;
	vector<pair<string, int>> name;
	Tree* child[70];
};

Tree * insert(Tree* tree, string word) {
	for (int i = 0; i < tree->name.size(); i++) 
		if (tree->name[i].first == word) 
			return tree->child[i];
	tree->name.push_back(make_pair(word, tree->name.size()));
	tree->child[tree->name.size() - 1] = new Tree(*tree);
	return tree->child[tree->name.size() - 1];

}

int f(pair<string, int>& a, pair<string, int>& b) {
	return a.first < b.first;
}

void show(Tree* tree) {
	for (int i = 0; i < tree->name.size(); i++) {
		sort(tree->name.begin(), tree->name.end(), f);
		for (int j = 0; j < tree->level; j++)cout << " ";
		cout << (tree->name[i].first) << endl;
		show(tree->child[tree->name[i].second]);
	}
}

int main() {
	int n, count = 0;
	cin >> n;
	string temp;
	Tree tree;
	Tree* tempTree;
	for (int i = 0; i < n; i++) {
		cin >> temp;
		string single = "";
		tempTree = &tree;
		for (int j = 0; j < temp.size(); j++) {
			if (temp[j] == '\\') {
				tempTree = insert(tempTree, single);
				single = "";
			}
			else single += temp[j];
		}
		tempTree = insert(tempTree, single);
	}
	show(&tree);
}

#include<iostream>

#include<algorithm>

#include<vector>

using namespace std;

class Tree {

public:

Tree() {

level = 0;

}

Tree(const Tree& tree) {

this->level = tree.level + 1;

}

int level;

vector<pair<string, int>> name;

Tree* child[70];

};

Tree * insert(Tree* tree, string word) {

for (int i = 0; i < tree->name.size(); i++)

if (tree->name[i].first == word)

return tree->child[i];

tree->name.push_back(make_pair(word, tree->name.size()));

tree->child[tree->name.size() - 1] = new Tree(*tree);

return tree->child[tree->name.size() - 1];

}

int f(pair<string, int>& a, pair<string, int>& b) {

return a.first < b.first;

}

void show(Tree* tree) {

for (int i = 0; i < tree->name.size(); i++) {

sort(tree->name.begin(), tree->name.end(), f);

for (int j = 0; j < tree->level; j++)cout << " ";

cout << (tree->name[i].first) << endl;

show(tree->child[tree->name[i].second]);

}

int main() {

int n, count = 0;

cin >> n;

string temp;

Tree tree;

Tree* tempTree;

for (int i = 0; i < n; i++) {

cin >> temp;

string single = "";

tempTree = &tree;

for (int j = 0; j < temp.size(); j++) {

if (temp[j] == '\\') {

tempTree = insert(tempTree, single);

single = "";

}

else single += temp[j];

}

tempTree = insert(tempTree, single);

}

show(&tree);

}