非线性方程求解

二分搜索法求解方程在给定区间内的实根

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
import math
import numpy as np
import random
def RootHalves(left,right,step,eps,result,m):
    n=0
    js=0
    z=left
    y=FunctionValue(z)
    while z<=(right+step/2) and n!=m:
        if abs(y)<eps:
            n=n+1
            result.append(z)
            z=z+step/2.0
            z=FunctionValue(z)
        else:
            z1=z+step
            y1=FunctionValue(z1)
            if(abs(y1)<eps):
                n=n+1
                result.append(z)
                z=z1+step/2
                y=FunctionValue(z)
            else:
                if y*y1>0.0:
                    y=y1
                    z=z1
                else:
                    js=0
                    while js==0:
                        if abs(z1-z)<eps:
                            n=n+1
                            result.append((z1+z)/2.0)
                            z=z1+step/2.0
                            y=FunctionValue(z)
                            js=1
                        else:
                            z0=(z1+z)/2.0
                            y0=FunctionValue(z0)
                            if abs(y0)<eps:
                                result[n]=0
                                n=n+1
                                js=1
                                z=z0+step/2.0
                                y=FunctionValue(z)
                            else:
                                if (y*y0)<0.0:
                                    z1=z0
                                    y1=y0
                                else:
                                    z=z0
                                    y=y0
    
    return(n)


对二分搜索法的一些测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
##求一维多项式值
def PolyValueOneDim(dCoff,stNo,dx):
    dValue=dCoff[stNo-1]
    for st in range(stNo-2,-1,-1):
        dValue=dValue*dx+dCoff[st]
    return(dValue)

def FunctionValue(dx):
    dCoff=[-20,7,-7,1,3,-5,1]
    dValue=PolyValueOneDim(dCoff,7,dx)
    return dValue


# 测试用例
left=-2
right=5
eps=1.e-8
step=0.2
result=[]
rootNum=RootHalves(left,right,step,eps,result,10)
print("rootNum = {}".format(rootNum))
print(result)
#rootNum = 2
#[-1.40246, 4.33375]
rootNum = 2
[-1.4024630278348929, 4.33375544846058]

牛顿（newton）法求解非线性方程的一个实根

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# x储存迭代初值，返回时存放迭代终值
# eps 控制精度
# js 最大迭代次数
# 返回值<0表示出现f'(x)=0，>=js表示精度不满足要求,>=0表示正常返回

def RootNewton(x,eps,js):
    k=1
    r=js
    x0=x
    y=FunctionValue(x0)
    y1=FunctionCValue(x0)
    d=eps+1.0
    while(d>eps or abs(d-eps)<eps)and (r!=0):
        if abs(y1)<eps:
            return(-1)
        else:
            x1=x0-y/y1
            y=FunctionValue(x1)
            y1=FunctionCValue(x1)
            d=abs(x1-x0)
            p=abs(y)
            if p>d:
                d=p
            x0=x1
            r-=1
    
    x=x1
    k=js-r
    return(k,x)

对牛顿法的一些测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#求函数的值
def FunctionValue(dx):
    dCoff=[-1,0,-1,1]
    dValue=PolyValueOneDim(dCoff,4,dx)
    return dValue

#求函数的导数值
def FunctionCValue(dx):
    dCoff1=[0,-2,3]
    dValue=PolyValueOneDim(dCoff1,3,dx)
    return dValue


x=1.5
print("------------RootNewton()---------")
k,x=RootNewton(x,1.e-6,60)
print("k={} x={}".format(k,x))
    
------------RootNewton()---------
k=60 x=1.465571231876768

埃特金(AitKen)法求解非线性方程的一个实根

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# x储存迭代初值，返回时存放迭代终值
# eps 控制精度
# js 最大迭代次数
# 返回值<0表示出现f'(x)=0，>=js表示精度不满足要求,>=0表示正常返

def RootAitken(x,eps,js):
    r=0
    flag=0
    x0=x
    while flag==0 and r!=js:
        r+=1
        u=FunctionValue(x0)
        v=FunctionValue(u)
        if abs(u-v)<eps:
            x0=v
            flag=1
        else:
            x0=v-(v-u)*(v-u)/(v-2.0*u+x0)
    
    x=x0
    r=js-r
    return (r,x)

对埃特金(AitKen)法的一些测试

1
2
3
4
5
6
7
8
9
#求函数的值
def FunctionValue(dx):
    dCoff=[6,0,-1]
    dValue=PolyValueOneDim(dCoff,3,dx)
    return dValue

r,x=RootAitken(0,1.e-6,20)
print("---------RootAitKen------------")
print("k={} x={}".format(r,x))
---------RootAitKen------------
k=13 x=1.9999994417336957

连分式(Fraction)法求解非线性方程的一个实根

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# x储存迭代初值，返回时存放迭代终值
# eps 控制精度
#返回值<10表示正常返回
def RootFraction(x,eps):
    r=10
    q=1.0e+35
    h=0
    a=[]
    y=[]
    x0=x
    while r!=0:
        r=r-1
        j=0
        it=r
        while j<=7:
            if j<=2:
                z=x0+0.1*j
            else:
                z=h
            y.append(FunctionValue(z))
            h=z
            if(j==0):
                a.append(z)
            else:
                m=0
                i=0
                while m==0 and i<j:
                    if abs(h-a[i])<=eps:
                        m=1
                    else:
                        h=(y[j]-y[i])/(h-a[i])
                    i+=1
                a.append(h)
                if m!=0:
                    a[j]=q
                h=0.0
                for i in range(j-1,-1,-1):
                    if abs(a[i+1]+h)<=eps:
                        h=q
                    else:
                        h-=y[i]/(a[i+1]+h)
                
                h=h+a[0]
            if abs(y[j])>=eps:
                j+=1
            else:
                j=10
                r=0
        x0=h
    x=h
    return(10-it,x)

对连分式(Fraction)法的一些测试

1
2
3
4
5
6
7
8
9
10
#求函数的值
def FunctionValue(dx):
    dCoff=[-1,0,-1,1]
    dValue=PolyValueOneDim(dCoff,4,dx)
    return dValue

x=1.0
print("-----------RootFraction()---------")
k,x=RootFraction(x,1.0e-6)
print("k={} x={}".format(k,x))
-----------RootFraction()---------
k=10 x=47.30994086939946

梯度(Gradient)法求解非线性方程组的一组实根

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# x 长度为n，存放一组初值x0,X1，X2.....，返回时存放一组实根
# eps 控制精度
# js 最大迭代次数

def RootGradient(eps,x,js):
    n=len(x)
    y=[x for x in range(n)]
    r=js
    f=FunctionValue(x,y)
    while f>=eps:
        r=r-1
        if r==0:
            return js
        d=0
        for j in range(n):
            d=d+y[j]*y[j]
        if abs(d)<=eps:
            return -1
        s=f/d
        for j in range(n):
            x[j]=x[j]-s*y[j]
        
        f=FunctionValue(x,y)
    return(js-r)

对梯度(Gradient)法的一些测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def FunctionValue(x,y):
    f1=x[0]-5.0*x[1]*x[1]+7.0*x[2]*x[2]+12.0 #计算二维多项式f1函数值
    f2=3.0*x[0]*x[1]+x[0]*x[2]-11.0*x[0] #计算二维多项式f2函数值
    f3=2.0*x[1]*x[2]+40.0*x[0] #计算二维多项式f3函数值
    z=f1*f1+f2*f2+f3*f3
    df1=1.0
    df2=3.0*x[1]+x[2]-11
    df3=40
    y[0]=2.0*(f1*df1+f2*df2+f3*df3)
    df1=10.0*x[1]
    df2=3.0*x[0]
    df3=2.0*x[2]
    y[1]=2.0*(f1*df1+f2*df2+f3*df3)
    df1=14.0*x[2]
    df2=x[0]
    df3=2.0*x[1];
    y[2]=2.0*(f1*df1+f2*df2+f3*df3)
    return z

js=500
x=[1.5,6.5,-5.0]
i=RootGradient(1.0e-6,x,js)
print("------------RootGradient()-------------")
if i>0 and i<js:
    for i in range(3):
        print("x({})={}".format(i,x[i]))

------------RootGradient()-------------
x(0)=1.0001884137008277
x(1)=5.00043177648948
x(2)=-4.00038848410507

蒙特卡洛(MonteCarlo)法求非线性方程的一个实数根

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
# x：存放初值，返回方程终值
# b：均匀分布随机数的端点初值
# m：控制调节点b的参数
# eps：控制精度要求

def RootMonteCarloReal(x,b,m,eps):
    a=b
    k=1
    r=1.0
    X=x
    y=FunctionValue(X)
    while a>eps or abs(a-eps)<=0:
        x1=random.uniform(1.0,5.0)
        x1=-a+2.0*a*x1
        x1=X+x1
        y1=FunctionValue(x1)
        k+=1
        if abs(y1)>abs(y) or abs(abs(y1)-abs(y))<=eps:
            if k>m:
                k=1
                a/=2.0
        else:
            k=1
            X=x1
            y=y1
            if abs(y)<eps:
                x=X
                break
    
    x=X
    return x

对 蒙特卡洛(MonteCarlo)法的一些测试

1
2
3
4
5
6
7
8
def FunctionValue(x):
    y=math.exp(-x*x*x)-math.sin(x)/math.cos(x)+800.00
    return y

x=0.5
x=RootMonteCarloReal(x,1,10,1.0e-6)
print("-----------RootMonteCarloReal()---------------")
print("x={}".format(x))
-----------RootMonteCarloReal()---------------
x=10.994575404027522

1