二分法解超越方程的应用

# -*- coding: utf-8 -*-
'''二分法'''

import numpy as np
d = 20  #弦长
s = 36.234  #面积
a=0.01; b=3.14;
global step
step=0

def f(x,d,s):
    return x-np.sin(x)-8*s/(d*d)*np.sin(x/2)*np.sin(x/2)

def erfen(a,b,d,s): #a、b为方程的根区间，d为弦长，s为面积
    global step
    fhalf=f((a+b)/2,d,s);
    half=(a+b)/2;
    fa=f(a,d,s);fb=f(b,d,s);
    step=step+1;
    if (fa==0 ):
        return a
    if (fb==0):
        return a
    if(fhalf==0):
        return fhalf
    if np.sqrt(abs(fa*fb))<1e-10:
        return a
    if fhalf*fb<0:
        return erfen(half,b,d,s)
    else:
        return erfen(a,half,d,s)

x=erfen(a,b,d,s) #x即为方程的解
r=np.sqrt(d*d/(1-np.cos(x))/2.0)
print(r,x)