Python – 第 2 页 – Pancake's Personal Website

Python 多线程与多进程

python多线程因为GIL锁的原因可以说是假的，但是多进程是可以使用多核 cpu。

测试代码如下：

from threading import Thread

def loop1():
    a = 0
    while 1:
        a += 1

def loop2(a):
    while 1:
        a += 1

for i in range(3):
    t = Thread(target=loop1())
    t.start()

for i in range(3):
    t = Thread(target=loop2(), args=(1,))
    t.start()

from threading import Thread

def loop1():

a = 0

while 1:

a += 1

def loop2(a):

while 1:

a += 1

for i in range(3):

t = Thread(target=loop1())

t.start()

for i in range(3):

t = Thread(target=loop2(), args=(1,))

t.start()

from multiprocessing import Process, Pool

def loop1():
    a = 0
    while 1:
        a += 1

def loop2(a):
    while 1:
        a += 1

pool = Pool(processes=6)

for i in range(3):
    pool.apply_async(loop1)

for i in range(3):
    pool.apply_async(loop2, args=(1,))

pool.close()
pool.join()

from multiprocessing import Process, Pool

def loop1():

a = 0

while 1:

a += 1

def loop2(a):

while 1:

a += 1

pool = Pool(processes=6)

for i in range(3):

pool.apply_async(loop1)

for i in range(3):

pool.apply_async(loop2, args=(1,))

pool.close()

pool.join()

测试结果如下：

结果：我们可以看到 threading 多线程库只可以跑满一核 cpu，但是 multiprocessing 多进程库可以跑满 cpu

二重积分参考代码

class doubleIntegral:
    def __init__(self, xMin, xMax, yMin, yMax, func):
        self.xMin = xMin
        self.xMax = xMax
        self.yMin = yMin
        self.yMax = yMax
        self.func = func

    def trapezoid(self, nX, nY):
        result = 0
        betweenX = (self.xMax - self.xMin) / nX
        betweenY = (self.yMax - self.yMin) / nY
        sums = []
        for i in range(nY + 1):
            temp = 0
            for j in range(nX):
                temp += betweenX / 2 * (self.func(self.xMin + j * betweenX, self.yMin + i * betweenY) +
                                        self.func(self.xMin + (j + 1) * betweenX, self.yMin + i * betweenY))
            sums.append(temp)
        ans = 0
        for i in range(nY):
            ans += betweenY / 2 * (sums[i] + sums[i + 1])
        return ans

    def simpsons(self, nX, nY):
        result = 0
        betweenX = (self.xMax - self.xMin) / nX
        betweenY = (self.yMax - self.yMin) / nY / 2
        sums = []
        for i in range(2 * nY + 1):
            temp = 0
            for j in range(nX):
                a = self.xMin + j * betweenX
                b = self.xMin + (j + 1) * betweenX
                c = (a + b) / 2
                temp += betweenX / 6 * (self.func(a, self.yMin + i * betweenY) + self.func(
                    b, self.yMin + i * betweenY) + 4 * self.func(c, self.yMin + i * betweenY))
            sums.append(temp)
        ans = 0
        for i in range(nY):
            ans += betweenY / 3 * \
                (sums[i * 2] + sums[(i + 1) * 2] + 4 * sums[i * 2 + 1])
        return ans


if __name__ == "__main__":
    # 个人不知道为什么梯形法求出来不对，有兴趣的可以自己去看看
    # 输入格式：xMin, xMax, yMin, yMax, func
    # 输入格式'： nx, ny
    print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).trapezoid(4, 8))
    print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).simpsons(4, 8))

class doubleIntegral:

def __init__(self, xMin, xMax, yMin, yMax, func):

self.xMin = xMin

self.xMax = xMax

self.yMin = yMin

self.yMax = yMax

self.func = func

def trapezoid(self, nX, nY):

result = 0

betweenX = (self.xMax - self.xMin) / nX

betweenY = (self.yMax - self.yMin) / nY

sums = []

for i in range(nY + 1):

temp = 0

for j in range(nX):

temp += betweenX / 2 * (self.func(self.xMin + j * betweenX, self.yMin + i * betweenY) +

self.func(self.xMin + (j + 1) * betweenX, self.yMin + i * betweenY))

sums.append(temp)

ans = 0

for i in range(nY):

ans += betweenY / 2 * (sums[i] + sums[i + 1])

return ans

def simpsons(self, nX, nY):

result = 0

betweenX = (self.xMax - self.xMin) / nX

betweenY = (self.yMax - self.yMin) / nY / 2

sums = []

for i in range(2 * nY + 1):

temp = 0

for j in range(nX):

a = self.xMin + j * betweenX

b = self.xMin + (j + 1) * betweenX

c = (a + b) / 2

temp += betweenX / 6 * (self.func(a, self.yMin + i * betweenY) + self.func(

b, self.yMin + i * betweenY) + 4 * self.func(c, self.yMin + i * betweenY))

sums.append(temp)

ans = 0

for i in range(nY):

ans += betweenY / 3 * \

(sums[i * 2] + sums[(i + 1) * 2] + 4 * sums[i * 2 + 1])

return ans

if __name__ == "__main__":

# 个人不知道为什么梯形法求出来不对，有兴趣的可以自己去看看

# 输入格式：xMin, xMax, yMin, yMax, func

# 输入格式'： nx, ny

print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).trapezoid(4, 8))

print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).simpsons(4, 8))

计算数学 ODE LAB 3

全部是正确答案，慎用。Adams和Milne需要用另外的算法如runge-kutta（代码示例求法）求前3个数，当然也可以用Euler和Euler Improved，大家不要都一样啦。Python缩进规范请注意！

Python Euler Improved

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

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

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

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

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

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

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

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

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]

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]

计算数学 ODE Python参考代码

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

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

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)