用 Python 解数学二分法题

瞎扯

不知各位上过高中没,我现在也才高一,数学学的有点痛苦,但是教科书里有不少关于信息技术实践的项目。虽然学校条件不足,老师不讲,但是有好奇心的我还是乐于去尝试的。

正文

今天就让大家见识一下 二分法(bisection method) 在 Python 中的实现方法。

流程图

这里直接放上《普通高中教科书 数学 必修 第一册(RJ)》的 P146 页,包含笔算和流程图,写的过程中也是需要配合看的,而且程序逻辑比笔算逻辑更严密。

据此,我们基本上可以开始用编程语言编写程序了。

一个小插叙

这是我在晚修时花了两分钟手写的草稿代码:

我拿出来展示只是为了让大家函数名要牢记,不要像我一样。请把目光从那个 绝对值() 上移开

导入库

考虑到上了高中后的数学对于对数(logarithm)有需要,我在代码里加了以下代码(注释掉了,如需使用请删去注释符号):

from math import log

# Usage/用法:
# log(antilogarithm, base)
# log(真数, 底数)
# Attention: The order of the antilogarithm and the base is different to the writing habits.
# 注意:真数和底数的顺序不同于笔算习惯。

定义数学函数

在第一个 def 代码块,我定义了一个名为 f 的函数,它有一个形参:x。在缩进的内容里,注释将引导你定义一个数学函数并使其可以返回函数值 y 的值。

定义的方式和数学很像,你只需要将数学方程式子转化为程序赋值式即可,当然最简单的就是没有 k 的正比例函数了,直接 y = x

def f(x):
    x = float(x)
    # You need to define the math function that you want to calculate here first.
    # The function expression should be a program expression, that is, shortcuts such as omitting the multiplication sign cannot be used.

    # y = 2**x + 3*x - 7  # exponential function test
    # y = log(x) + 2 * x - 6  # logarithmic function test
    # y = x ** 3 + 1.1 * x ** 2 + 0.9 * x - 1.4  # A question in textbook..
    # y = 1 * x  # linear function test
    # y = 3 - log(x, 10) - x  # A question in textbook
    # y = log(x) - 2 / x  # A question in textbook

    # Example:
    y = x

    return y

如果你不会将数学表达式转为程序赋值式,我只能认为你没上过学,至少你的学校经费不足。

主函数

在主函数里,我根据流程图的内容先设置了三个输入,分别获取精确度(ε),区间头(a),区间尾(b)。

还做了取模处理,让等于整数的浮点数直接显示为整数。

epsilon = float(input("ε= "))  # Accuracy
a = float(input("a="))  # Start of range
b = float(input("b="))  # End of range

if a % 1 == 0:
    a = int(a)
if b % 1 == 0:
    b = int(b)

对于中间的计算和判断,我用了一个 while 循环,如果条件相反则一直循环,直到条件不符那个相反的条件就跳出循环输出答案。其实说简单点就是直接把流程图翻译为代码。

while True:
    c = (a+b) / 2
    if c % 1 == 0:
        c = int(c)
    else:
        pass

    if f(a) * f(c):
        b = c
    else:
        if f(c) == 0:
            a = c
            break
        else:
            a = c

跳出循环后,就只有输出结果了,再加上一点人性化处理,就是下面这样:

x = a
if x % 1 == 0:
    x = int(x)
else:
    pass
print("x=" + str(x))

效果

程序写到这,基本上运行起来就已经很正常了,就是在你输完三个参数后,回车,结果就出来了(速度取决于电脑计算速度和精确度的大小)

# y = x - 1

[Input]
ε= 1
a=0
b=2

[Output]
x=1

改进

对于此程序,用来为了直接获得结果的计算已经够用,但是中途我为了调试和笔算验证,便使用了一些输出,并让输出呈现为数学证明过程,效果大致如下:

这里拿的是课本里的练习题,最终用作业帮查答案验证:

例1

# y = x ** 3 + 1.1 * x ** 2 + 0.9 * x - 1.4

[Input]
ε= 0.1
a=0
b=1

[Output]
c=(0+1)/2=0.5
∵f(b)f(c)=1.6×-0.5499999999999998=-0.8799999999999998<0
∴a=c=0.5
∴x₀∈(0.5,1)
c=(0.5+1)/2=0.75
∵f(a)f(c)=-0.5499999999999998×0.31562500000000004=-0.17359374999999996<0
∴b=c=0.75
∴x₀∈(0.5,0.75)
c=(0.5+0.75)/2=0.625
∵f(b)f(c)=0.31562500000000004×-0.1636718749999999=-0.05165893554687498<0
∴a=c=0.625
∴x₀∈(0.625,0.75)
c=(0.625+0.75)/2=0.6875
∵f(a)f(c)=-0.1636718749999999×0.06362304687500009=-0.010413303375244149<0
∴b=c=0.6875
∴x₀∈(0.625,0.6875)
∵|0.625-0.6875|<0.1
∴x=a=0.625| < 0.1
∴x=a=0.625

由于遇到了一些二进制中的无限小数,所以在计算过程中输出了一个不太对劲的数,但是比较逼近精确值,所以最终还是把答案算出来了。

验证,答案正确。

但是作为一个辅助运算的东西,这个东西必须抗压,所以当我把精确度设置为 0.0000000000000000001 后,这个东西也是成功帮我把结果算到了一个更加精确的数字:

例2

[Output]
c=(0+1)/2=0.5
∵f(b)f(c)=1.6×-0.5499999999999998=-0.8799999999999998<0
∴a=c=0.5
∴x₀∈(0.5,1)
c=(0.5+1)/2=0.75
∵f(a)f(c)=-0.5499999999999998×0.31562500000000004=-0.17359374999999996<0
∴b=c=0.75
∴x₀∈(0.5,0.75)
c=(0.5+0.75)/2=0.625
∵f(b)f(c)=0.31562500000000004×-0.1636718749999999=-0.05165893554687498<0
∴a=c=0.625
∴x₀∈(0.625,0.75)
c=(0.625+0.75)/2=0.6875
∵f(a)f(c)=-0.1636718749999999×0.06362304687500009=-0.010413303375244149<0
∴b=c=0.6875
∴x₀∈(0.625,0.6875)
c=(0.625+0.6875)/2=0.65625
∵f(b)f(c)=0.06362304687500009×-0.05302124023437482=-0.00337337285280227<0
∴a=c=0.65625
∴x₀∈(0.65625,0.6875)
c=(0.65625+0.6875)/2=0.671875
∵f(a)f(c)=-0.05302124023437482×0.0045402526855471415=-0.00024072982836516045<0
∴b=c=0.671875
∴x₀∈(0.65625,0.671875)
c=(0.65625+0.671875)/2=0.6640625
∵f(b)f(c)=0.0045402526855471415×-0.02442922592163077=-0.00011091485859652195<0
∴a=c=0.6640625
∴x₀∈(0.6640625,0.671875)
c=(0.6640625+0.671875)/2=0.66796875
∵f(b)f(c)=0.0045402526855471415×-0.009991848468780429=-4.5365516843960435e-05<0
∴a=c=0.66796875
∴x₀∈(0.66796875,0.671875)
c=(0.66796875+0.671875)/2=0.669921875
∵f(b)f(c)=0.0045402526855471415×-0.002737660706042977=-1.242967137272851e-05<0
∴a=c=0.669921875
∴x₀∈(0.669921875,0.671875)
c=(0.669921875+0.671875)/2=0.6708984375
∵f(a)f(c)=-0.002737660706042977×0.0008983274921776641=-2.4593158764929206e-06<0
∴b=c=0.6708984375
∴x₀∈(0.669921875,0.6708984375)
c=(0.669921875+0.6708984375)/2=0.67041015625
∵f(b)f(c)=0.0008983274921776641×-0.0009204083820804065=-8.268281536535929e-07<0
∴a=c=0.67041015625
∴x₀∈(0.67041015625,0.6708984375)
c=(0.67041015625+0.6708984375)/2=0.670654296875
∵f(b)f(c)=0.0008983274921776641×-1.1225932393887916e-05=-1.0084563694757333e-08<0
∴a=c=0.670654296875
∴x₀∈(0.670654296875,0.6708984375)
c=(0.670654296875+0.6708984375)/2=0.6707763671875
∵f(a)f(c)=-1.1225932393887916e-05×0.00044350440257412416=-4.978750439688768e-09<0
∴b=c=0.6707763671875
∴x₀∈(0.670654296875,0.6707763671875)
c=(0.670654296875+0.6707763671875)/2=0.67071533203125
∵f(a)f(c)=-1.1225932393887916e-05×0.00021612764144274266=-2.4262342912866775e-09<0
∴b=c=0.67071533203125
∴x₀∈(0.670654296875,0.67071533203125)
c=(0.670654296875+0.67071533203125)/2=0.670684814453125
∵f(a)f(c)=-1.1225932393887916e-05×0.00010244795619795966=-1.1500738301702857e-09<0
∴b=c=0.670684814453125
∴x₀∈(0.670654296875,0.670684814453125)
c=(0.670654296875+0.670684814453125)/2=0.6706695556640625
∵f(a)f(c)=-1.1225932393887916e-05×4.561028733096606e-05=-5.120180020432276e-10<0
∴b=c=0.6706695556640625
∴x₀∈(0.670654296875,0.6706695556640625)
c=(0.670654296875+0.6706695556640625)/2=0.6706619262695312
∵f(a)f(c)=-1.1225932393887916e-05×1.7191996326992864e-05=-1.9299618848279128e-10<0
∴b=c=0.6706619262695312
∴x₀∈(0.670654296875,0.6706619262695312)
c=(0.670654296875+0.6706619262695312)/2=0.6706581115722656
∵f(a)f(c)=-1.1225932393887916e-05×2.9829866812214334e-06=-3.3486806815259896e-11<0
∴b=c=0.6706581115722656
∴x₀∈(0.670654296875,0.6706581115722656)
c=(0.670654296875+0.6706581115722656)/2=0.6706562042236328
∵f(b)f(c)=2.9829866812214334e-06×-4.121484177499468e-06=-1.2294332408345787e-11<0
∴a=c=0.6706562042236328
∴x₀∈(0.6706562042236328,0.6706581115722656)
c=(0.6706562042236328+0.6706581115722656)/2=0.6706571578979492
∵f(b)f(c)=2.9829866812214334e-06×-5.692515785415964e-07=-1.6980698770538587e-12<0
∴a=c=0.6706571578979492
∴x₀∈(0.6706571578979492,0.6706581115722656)
c=(0.6706571578979492+0.6706581115722656)/2=0.6706576347351074
∵f(a)f(c)=-5.692515785415964e-07×1.206866843794785e-06=-6.870108559196955e-13<0
∴b=c=0.6706576347351074
∴x₀∈(0.6706571578979492,0.6706576347351074)
c=(0.6706571578979492+0.6706576347351074)/2=0.6706573963165283
∵f(a)f(c)=-5.692515785415964e-07×3.1880745576806646e-07=-1.81481647446802e-13<0
∴b=c=0.6706573963165283
∴x₀∈(0.6706571578979492,0.6706573963165283)
c=(0.6706571578979492+0.6706573963165283)/2=0.6706572771072388
∵f(b)f(c)=3.1880745576806646e-07×-1.2522210557364133e-07=-3.992174088385281e-14<0
∴a=c=0.6706572771072388
∴x₀∈(0.6706572771072388,0.6706573963165283)
c=(0.6706572771072388+0.6706573963165283)/2=0.6706573367118835
∵f(a)f(c)=-1.2522210557364133e-07×9.679266410600462e-08=-1.2120581203436115e-14<0
∴b=c=0.6706573367118835
∴x₀∈(0.6706572771072388,0.6706573367118835)
c=(0.6706572771072388+0.6706573367118835)/2=0.6706573069095612
∵f(b)f(c)=9.679266410600462e-08×-1.4214723620398217e-08=-1.3758809687488946e-15<0
∴a=c=0.6706573069095612
∴x₀∈(0.6706573069095612,0.6706573367118835)
c=(0.6706573069095612+0.6706573367118835)/2=0.6706573218107224
∵f(a)f(c)=-1.4214723620398217e-08×4.128896957666939e-08=-5.869112911033857e-16<0
∴b=c=0.6706573218107224
∴x₀∈(0.6706573069095612,0.6706573218107224)
c=(0.6706573069095612+0.6706573218107224)/2=0.6706573143601418
∵f(a)f(c)=-1.4214723620398217e-08×1.3537122978135585e-08=-1.9242646174953936e-16<0
∴b=c=0.6706573143601418
∴x₀∈(0.6706573069095612,0.6706573143601418)
c=(0.6706573069095612+0.6706573143601418)/2=0.6706573106348515
∵f(b)f(c)=1.3537122978135585e-08×-3.388005431759211e-10=-4.586384618031579e-18<0
∴a=c=0.6706573106348515
∴x₀∈(0.6706573106348515,0.6706573143601418)
c=(0.6706573106348515+0.6706573143601418)/2=0.6706573124974966
∵f(a)f(c)=-3.388005431759211e-10×6.599161217479832e-09=-2.23579940498764e-18<0
∴b=c=0.6706573124974966
∴x₀∈(0.6706573106348515,0.6706573124974966)
c=(0.6706573106348515+0.6706573124974966)/2=0.670657311566174
∵f(a)f(c)=-3.388005431759211e-10×3.1301803371519554e-09=-1.0605067984656704e-18<0
∴b=c=0.670657311566174
∴x₀∈(0.6706573106348515,0.670657311566174)
c=(0.6706573106348515+0.670657311566174)/2=0.6706573111005127
∵f(a)f(c)=-3.388005431759211e-10×1.3956902300549245e-09=-4.728606080479347e-19<0
∴b=c=0.6706573111005127
∴x₀∈(0.6706573106348515,0.6706573111005127)
c=(0.6706573106348515+0.6706573111005127)/2=0.6706573108676821
∵f(a)f(c)=-3.388005431759211e-10×5.284448434395017e-10=-1.7903739999581778e-19<0
∴b=c=0.6706573108676821
∴x₀∈(0.6706573106348515,0.6706573108676821)
c=(0.6706573106348515+0.6706573108676821)/2=0.6706573107512668
∵f(a)f(c)=-3.388005431759211e-10×9.482237217639522e-11=-3.212587119859205e-20<0
∴b=c=0.6706573107512668
∴x₀∈(0.6706573106348515,0.6706573107512668)
c=(0.6706573106348515+0.6706573107512668)/2=0.6706573106930591
∵f(b)f(c)=9.482237217639522e-11×-1.2198908549976295e-10=-1.156729446671662e-20<0
∴a=c=0.6706573106930591
∴x₀∈(0.6706573106930591,0.6706573107512668)
c=(0.6706573106930591+0.6706573107512668)/2=0.670657310722163
∵f(b)f(c)=9.482237217639522e-11×-1.3583356661683865e-11=-1.2880061007789048e-21<0
∴a=c=0.670657310722163
∴x₀∈(0.670657310722163,0.6706573107512668)
c=(0.670657310722163+0.6706573107512668)/2=0.6706573107367149
∵f(a)f(c)=-1.3583356661683865e-11×4.061950775735568e-11=-5.517492612901967e-22<0
∴b=c=0.6706573107367149
∴x₀∈(0.670657310722163,0.6706573107367149)
c=(0.670657310722163+0.6706573107367149)/2=0.6706573107294389
∵f(a)f(c)=-1.3583356661683865e-11×1.3518075547835906e-11=-1.8362084154584262e-22<0
∴b=c=0.6706573107294389
∴x₀∈(0.670657310722163,0.6706573107294389)
c=(0.670657310722163+0.6706573107294389)/2=0.6706573107258009
∵f(b)f(c)=1.3518075547835906e-11×-3.2862601528904634e-14=-4.442391301661606e-25<0
∴a=c=0.6706573107258009
∴x₀∈(0.6706573107258009,0.6706573107294389)
c=(0.6706573107258009+0.6706573107294389)/2=0.6706573107276199
∵f(a)f(c)=-3.2862601528904634e-14×6.742606473153501e-12=-2.215795897934565e-25<0
∴b=c=0.6706573107276199
∴x₀∈(0.6706573107258009,0.6706573107276199)
c=(0.6706573107258009+0.6706573107276199)/2=0.6706573107267104
∵f(a)f(c)=-3.2862601528904634e-14×3.355093980417223e-12=-1.1025711657047777e-25<0
∴b=c=0.6706573107267104
∴x₀∈(0.6706573107258009,0.6706573107267104)
c=(0.6706573107258009+0.6706573107267104)/2=0.6706573107262557
∵f(a)f(c)=-3.2862601528904634e-14×1.6611156894441592e-12=-5.45885829956151e-26<0
∴b=c=0.6706573107262557
∴x₀∈(0.6706573107258009,0.6706573107262557)
c=(0.6706573107258009+0.6706573107262557)/2=0.6706573107260283
∵f(a)f(c)=-3.2862601528904634e-14×8.142375662600898e-13=-2.6757964689870415e-26<0
∴b=c=0.6706573107260283
∴x₀∈(0.6706573107258009,0.6706573107260283)
c=(0.6706573107258009+0.6706573107260283)/2=0.6706573107259146
∵f(a)f(c)=-3.2862601528904634e-14×3.907985046680551e-13=-1.2842655536998072e-26<0
∴b=c=0.6706573107259146
∴x₀∈(0.6706573107258009,0.6706573107259146)
c=(0.6706573107258009+0.6706573107259146)/2=0.6706573107258578
∵f(a)f(c)=-3.2862601528904634e-14×1.7896795156957523e-13=-5.8813524788752535e-27<0
∴b=c=0.6706573107258578
∴x₀∈(0.6706573107258009,0.6706573107258578)
c=(0.6706573107258009+0.6706573107258578)/2=0.6706573107258293
∵f(a)f(c)=-3.2862601528904634e-14×7.327471962526033e-14=-2.4079979131871385e-27<0
∴b=c=0.6706573107258293
∴x₀∈(0.6706573107258009,0.6706573107258293)
c=(0.6706573107258009+0.6706573107258293)/2=0.6706573107258151
∵f(a)f(c)=-3.2862601528904634e-14×2.020605904817785e-14=-6.640236669697867e-28<0
∴b=c=0.6706573107258151
∴x₀∈(0.6706573107258009,0.6706573107258151)
c=(0.6706573107258009+0.6706573107258151)/2=0.670657310725808
∵f(b)f(c)=2.020605904817785e-14×-6.217248937900877e-15=-1.2562609915644613e-28<0
∴a=c=0.670657310725808
∴x₀∈(0.670657310725808,0.6706573107258151)
c=(0.670657310725808+0.6706573107258151)/2=0.6706573107258116
∵f(a)f(c)=-6.217248937900877e-15×7.105427357601002e-15=-4.417621069237666e-29<0
∴b=c=0.6706573107258116
∴x₀∈(0.670657310725808,0.6706573107258116)
c=(0.670657310725808+0.6706573107258116)/2=0.6706573107258098
∵f(a)f(c)=-6.217248937900877e-15×4.440892098500626e-16=-2.7610131682735413e-30<0
∴b=c=0.6706573107258098
∴x₀∈(0.670657310725808,0.6706573107258098)
c=(0.670657310725808+0.6706573107258098)/2=0.6706573107258089
∵f(b)f(c)=4.440892098500626e-16×-2.886579864025407e-15=-1.2818989709841442e-30<0
∴a=c=0.6706573107258089
∴x₀∈(0.6706573107258089,0.6706573107258098)
c=(0.6706573107258089+0.6706573107258098)/2=0.6706573107258094
∵f(b)f(c)=4.440892098500626e-16×-1.3322676295501878e-15=-5.9164567891575885e-31<0
∴a=c=0.6706573107258094
∴x₀∈(0.6706573107258094,0.6706573107258098)
c=(0.6706573107258094+0.6706573107258098)/2=0.6706573107258096
∵f(b)f(c)=4.440892098500626e-16×-4.440892098500626e-16=-1.9721522630525295e-31<0
∴a=c=0.6706573107258096
∴x₀∈(0.6706573107258096,0.6706573107258098)
c=(0.6706573107258096+0.6706573107258098)/2=0.6706573107258097
∵f(c)=0
∴a=c=0.6706573107258097
∵|0.6706573107258097-0.6706573107258098|<1e-19
∴x=a=0.6706573107258097

而这一切,只用了 2.5s 。(我的电脑为 华硕天选2

总结

很好,充分发挥了我的添砖加瓦能力。至于源码,篇幅有限,还是不展示了。如果你只是需要知道结果,那么你只需要将上面所出现的代码块自定义组装,运行即可。

直接让两位大哥帮忙托管:

Math/bisection.py at main · gkcoll/Math · GitHub

bisection.py · 灰尘极客/Math - Gitee.com

再来一个网盘备份:

bisection.py


补充说明

二分法求零点是在难以直接用 y=0 算出结果时才使用的,且要求图像在 y=0 附近连续不断;再这种方法算出来的只是估值,除非你把精确度设置得非常小(难以想象的那种小);最后你得把范围 ab 输对,否则会解除错误答案。

综上,

故建议大家能算出精确值时不要过度依赖二分法。

如何使用 Python 解方程(组)?

你可以试试使用 SumPy 库,参考链接

版权声明:
作者:灰尘疾客
链接:https://www.gkcoll.xyz/418.html
来源:极客藏源
文章版权归作者所有,未经允许请勿转载。

THE END
分享
二维码
< <上一篇
下一篇>>