from matplotlib.pyplot import *
from numpy import *

## Exercice 1

f1=lambda x:x**2-2
f2=lambda x:x**7-1000
f3=lambda x:4*x**4-2*x**3+x**2-18*x+8

def dichotomie(f,a,b,eps=2**(-51)):
    x=a
    y=b
    n=0
    while abs(x-y)>eps:
        n+=1
        z=(x+y)/2
        if f(x)*f(z)>0:
            x=z
        else:
            y=z
    return (x,n)

# print(dichotomie(f1,1,2,2**(-51)))
# print(dichotomie(f2,1,10,0.001))
# print(dichotomie(f3,0,1,0.001))

## Exercice 2

def newton(f,df,u0,eps=2**(-51)):
    h=lambda x:x-f(x)/df(x)
    u=u0
    v=h(u)
    n=0
    while abs(u-v)>eps:
        n+=1
        u=v
        v=h(v)
    return(u,n)

df1=lambda x:2*x
df2=lambda x:7*x**6
df3=lambda x:16*x**3-6*x**2+2*x-18

# print(newton(f1,df1,3,2**(-51)))
# print(newton(f2,df2,3,0.001))
# print(newton(f3,df3,0.5,0.001))

## Exercice 3

def Der(f,x,h=2**(-17)):
    return (f(x+h)-f(x-h))/(2*h)

## Exercice 4

def secante(f,u0,u1,eps=2**(-51)):
    u=u0
    v=u1
    n=0
    while abs(u-v)>eps:
        n+=1
        a=(f(v)-f(u))/(v-u)
        u=v
        v=v-f(v)/a
    return(u,n)

# print(secante(f1,0,3))
# print(secante(f2,0,3))
# print(secante(f3,0,3))

# ----

c=1
t=0.17

fa=lambda x:t*c*(1.4845*sqrt(x/c)-0.63*x/c-1.758*(x/c)**2+1.4215*(x/c)**3-0.5075*(x/c)**4)

dfa=lambda x:Der(fa,x)

# print(dichotomie(dfa,0.1,1))

## ----

a=-41.9
b=-359.4
c=15211
d=141150
fc=lambda p:p**4+a*p**3+b*p**2+c*p+d

# Bornes pour x :
x0,x1=-20,50
# Bornes pour y :
y0,y1=-50,50

x=linspace(x0,x1,500)
plot(x,fc(x))

# Axes :
grid()
axhline(color="black")
axvline(color="black")

# axis([x0,x1,y0,y1])

show()

# print(dichotomie(fc,-100,100))

## ----

def bassins(z):
    f=lambda x:x**3-1
    df=lambda x:3*x**2
    l=newton(f,df,z)[0]
    c=l.imag
    if abs(c)<0.01:
        return 0
    elif abs(c-sqrt(3)/2)<0.01:
        return 1
    else:
        return 2

def bassins2(z):
    f=lambda x:x**6-1
    df=lambda x:6*x**5
    l=newton(f,df,z)[0]
    c=l.imag
    d=l.real
    if abs(c)<0.01:
        if abs(d-1)<0.01:
            return 0
        else:
            return 1
    elif abs(c-sqrt(3)/2)<0.01:
        if abs(d-1/2)<0.01:
            return 2
        else:
            return 3
    else:
        if abs(d-1/2)<0.01:
            return 4
        else:
            return 5

print(bassins(complex(2,4)))

A=array([[0,1,2],[2,1,1]])
imshow(A,interpolation="nearest")
B=[arange(10)]
imshow(B,interpolation="nearest")

a,b,c,d=-1.5,1.5,-1.5,1.5
N=1000
X=linspace(a,b,N)
Y=linspace(c,d,N)
M=zeros([N,N])
for i in range(N):
    for j in range(N):
        M[i,j]=bassins2(complex(X[i],Y[j]))
imshow(M,interpolation="nearest")
show()