python求方程近似解_二分法解非线性方程的算法及Python代码

算法：

例1：寻找 $$f(x) = x^3 - 11.0 x^2 + 38.8 x - 41.77 $$的有根区间.# -*- coding: utf-8 -*-

"""

寻找f(x) = x**3 - 11.0 * x**2 + 38.8*x - 41.77的有根区间

@author: morxio

"""

import numpy as np

def fun1(x):

#方程等号左边函数f(x)

return x**3 - 11.0 * x**2 + 38.8*x - 41.77

def sign(val):

#符号函数

if val < 0: s = "-"

elif val > 0: s = "+"

else: s = "0"

return s

def main():

#打印f(x)的符号

for x in range(0, 7, 1):

val = fun1(x)

print( "%3d %8.3f %5s" % (x, val, sign(val)) )

if __name__ == '__main__':

main()

结果：0 -41.770 -

1 -12.970 -

2 -0.170 -

3 2.630 +

4 1.430 +

5 2.230 +

6 11.030 +

例2：使用二分法求方程$$f(x) = x^3 - x -1 = 0$$在区间$$[1,2]$$内的实根，要求精确到$$10^{-2}$$.# -*- coding: utf-8 -*-

"""

使用二分法求方程f(x) = x**3 - x -1 = 0在区间[1,2]内的实根，要求精确到10**-2

@author: morxio

"""

import numpy as np

def fun2(x):

return (x**3 - x -1)

def midVal(x1, x2):

return ( 0.5 * (x1 + x2) )

#二分法程序

def binarySolver(a0, b0, EPS, ETA):

a = np.array([a0])

b = np.array([b0])

x = np.array([midVal(a,b)])

#判断x0是否是解

if np.abs( fun2(x[0]) ) < ETA:

return (x[0], 0)

N = 100 #最大迭代步

for k in range(0,N,1):

print("第%d步: 近似解%.4f在区间[%.4f, %.4f]内, 函数值为%.4f"\

% (k, x[k], a[k], b[k], fun2(x[k])) )

fL = fun2(a[k])

fR = fun2(b[k])

fx = fun2(x[k])

#决定有根子区间

if fL * fx < 0:

a = np.append(a, a[k])

b = np.append(b, x[k])

x = np.append(x, midVal(a[k+1], b[k+1]))

elif (fx * fR < 0):

a = np.append(a, x[k])

b = np.append(b, b[k])

x = np.append(x, midVal(a[k+1], b[k+1]))

if np.abs(fun2(x[k+1])) < ETA:

return (x[k+1], k+1)

if np.abs(fun2(a[k+1])) < ETA:

return (a[k+1], k)

if np.abs(fun2(b[k+1])) < ETA:

return (b[k+1], k)

if np.abs(a[k+1] - b[k]) < EPS:

print("%d: %.4f in [%.4f, %.4f], 函数值为%.4f"\

% (k+1, x[k+1], a[k+1], b[k+1], fun2(x[k+1])) )

print("达到停止准则: b[%d] - a[%d] = %.4f < EPS"\

%(k+1, k+1, b[k+1] - a[k+1]) )

return (x[k+1], k+1)

#主函数

def main():

EPS = 1e-2 #近似解的误差限

ETA = 1e-6 #函数值误差限

a, b = 1., 2.

fa, fb = fun2(a), fun2(b)

if np.abs(fa) < ETA:

print(f"f(1) = {fa:.8f} => solution at x = 1")

elif np.abs(fb) < ETA:

print(f"f(2) = {fb:.8f} => solution at x = 2")

else:

#手动判断一阶导数大于零，因此有唯一解

(x,n) = binarySolver(a, b, EPS, ETA)

print(f"{n:d}次迭代后, 找到近似解: {x:.4f}，此时函数值为{fun2(x):.4f}")

if __name__ == '__main__':

main()

运行结果：第0步: 近似解1.5000在区间[1.0000, 2.0000]内, 函数值为0.8750

第1步: 近似解1.2500在区间[1.0000, 1.5000]内, 函数值为-0.2969

第2步: 近似解1.3750在区间[1.2500, 1.5000]内, 函数值为0.2246

第3步: 近似解1.3125在区间[1.2500, 1.3750]内, 函数值为-0.0515

第4步: 近似解1.3438在区间[1.3125, 1.3750]内, 函数值为0.0826

第5步: 近似解1.3281在区间[1.3125, 1.3438]内, 函数值为0.0146

第6步: 近似解1.3203在区间[1.3125, 1.3281]内, 函数值为-0.0187

7: 1.3242 in [1.3203, 1.3281], 函数值为-0.0021

达到停止准则: b[7] - a[7] = 0.0078 < EPS

7次迭代后, 找到近似解: 1.3242，此时函数值为-0.0021