science-computation/Lab7_newtons.py at master · raymondWai/science-computation · GitHub

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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import numpy as np
import matplotlib.pylab as plt
import Lab6_nonlinear_relaxation as relax

def binary_search(f,interval=[0,10],i=0,epi=10e-7):
    x=float(interval[0])+interval[1]
    x/=2
    df=f(x+epi)-f(x-epi)
    df/=epi
    if f(x)+epi<0:
        if df>0:
            return binary_search(f,[x,interval[1]],i+1,epi)
        elif df<0:
            return binary_search(f,[interval[0],x],i+1,epi)
        else:
            df=f(x+2*epi)-f(x)
            df/=epi
            if df>0:
                return binary_search(f,[x,interval[1]],i+1,epi)
            elif df<0:
                return binary_search(f,[interval[0],x],i+1,epi)
            else:
                return False
    elif f(x)-epi>0:
            if df<0:
                return binary_search(f,[x,interval[1]],i+1,epi)
            elif df>0:
                return binary_search(f,[interval[0],x],i+1,epi)
            else:
                df=f(x+2*epi)-f(x)
                df/=epi
                if df<0:
                    return binary_search(f,[x,interval[1]],i+1,epi)
                elif df>0:
                    return binary_search(f,[interval[0],x],i+1,epi)
                else:
                    return False
    elif np.abs(f(x))<epi:
        return x,i
def newton(f,epi,degree,interval):
    x=interval[0]
    df=f(x+1e-5)-f(x-1e-5)
    df/=(2e-5)
    d2f=f(x+1e-5)-2*f(x)+f(x-1e-5)
    d2f/=1e-10
    tmp=[x,x-f(x)/df]
    err=np.abs(tmp[1]-tmp[0])
    step=[]
    ans=[]
    no_of_root=0
    while no_of_root<degree :
        step.append(0)
        step[-1]+=1
        err=np.abs(tmp[1]-tmp[0])
        print(err)
        while err>epi:
            tmp[0]=tmp[1]
            df=f(tmp[0]+1e-5)-f(tmp[0]-1e-5)
            df/=(2e-5)
            d2f=f(tmp[0]+epi)-2*f(tmp[0])+f(tmp[0]-epi)
            d2f/=1e-10
            try:
                print(df,d2f)
            except ValueError:
                return False
            tmp[1]=tmp[0]-f(tmp[0])/df
            step[-1]+=1
            err=np.abs(((-1.*d2f)/(2.*df)))*err**2
            if np.abs(tmp[1])>1e200:
                ans.append("NA")
        if len(ans)>no_of_root:
            pass
        else:
            ans.append(tmp[1])
        no_of_root+=1
        x+=(interval[0]+interval[1])/float(degree)
        tmp[0]=x
        df=f(tmp[0]+epi)-f(tmp[0]-epi)
        df/=(2*epi)
        if df!=0:
            tmp[1]=tmp[0]-f(tmp[0])/df
        else:
            x+=epi
            tmp[0]=x
            df=f(tmp[0]+epi)-f(tmp[0]-epi)
            df/=epi
            tmp[1]=tmp[0]-f(tmp[0])/df
    return ans,step
def f_1(x):
    return 5*np.math.exp(-1*x)+x-5
def f_1_relax(x):
    return 5-5*np.exp(-1*x)
def polynomial(x):
    return 924*x**6-2772*x**5+3150*x**4-1680*x**3+420*x**2-42*x+1
'''
root,i=binary_search(f_1,[0,10])
print("Binary search method: ",root,i)
root,i=relax.relaxation(f_1_relax,10,10e-7)
print("Relaxation",root,i)
ans,step=newton(f_1,10e-7,2,[0,10])
print("newton's method:",ans,step)
x_axis=np.linspace(0,1)
y=[polynomial(x) for x in x_axis]
plt.plot(x_axis,y)
plt.show()
ans,step=newton(polynomial,10e-12,6,[0,1])
print("newton's method:",ans,step)
ans,step=binary_search(polynomial,[0,1],0,10e-12)
print("Binary search:",ans,step)
'''