arc_length.py

import numpy as np
import matplotlib.pyplot as plt

#usei _ antes do nome para indicar que sao "variaveis globais"
_inicio, _meio, _fim = 1.0, (1.0+5.0)/2.0, 5.0


def F(x, y, r=2.0):
    return x - 0.5 * np.cos(x) - r * np.arccos(1 - y / r) + np.sqrt(y * (2 * r - y))


def dF(x, y, r=2.0):
    return (r - y) / np.sqrt(y * (2.0 * r - y)) - 1 / np.sqrt(y * (2.0 * r - y)/r**2)

def newton(y_inicial=None, n_iter=5, dominio_x=(_inicio, _fim), n_itens=501):
    valores_x = np.linspace(*(dominio_x), num=n_itens)
    if y_inicial is None:
        valores_y = np.full(n_itens, 2.0)
    else:
        valores_y = np.array([y_inicial])
    for i, x in enumerate(valores_x):
        if(i > 0):
            valores_y[i] = valores_y[i-1]
        for j in range(n_iter):
            valores_y[i] -= F(x, valores_y[i])/dF(x, valores_y[i])
    return valores_x, valores_y


#derivada de três pontos centrada
def df_to_determine_suitable_h_centered( x ):
    vet_dhf = np.zeros((10,))

    for i in range(10):
        y_ini = valores_f[x]

        x0 = x - vet_h[i]
        x1 = x
        x2 = x + vet_h[i]

        x0, y0 = newton(y_inicial=y_ini, dominio_x=(x0, x0), n_itens=1)
        x1, y1 = newton(y_inicial=y_ini, dominio_x=(x1, x1), n_itens=1)
        x2, y2 = newton(y_inicial=y_ini, dominio_x=(x2, x2), n_itens=1)

        Dhf = (y2 - y0) / (2 * vet_h[i])
        vet_dhf[i] = Dhf

    return vet_dhf

#derivada de três pontos avançada
def df_to_determine_suitable_h_forward( x ):
    vet_dhf = np.zeros((10,))

    for i in range(10):
        y_ini = valores_f[x]

        x0 = x
        x1 = x + vet_h[i]
        x2 = x + 2*vet_h[i]

        x0, y0 = newton(y_inicial=y_ini, dominio_x=(x0, x0), n_itens=1)
        x1, y1 = newton(y_inicial=y_ini, dominio_x=(x1, x1), n_itens=1)
        x2, y2 = newton(y_inicial=y_ini, dominio_x=(x2, x2), n_itens=1)

        Dhf = ( -3*y0 + 4*y1 - y2 ) / (2 * vet_h[i])
        vet_dhf[i] = Dhf

    return vet_dhf

#derivada de três pontos atrasada
def df_to_determine_suitable_h_backward( x ):
    vet_dhf = np.zeros((10,))

    for i in range(10):
        y_ini = valores_f[x]

        x0 = x - 2*vet_h[i]
        x1 = x - vet_h[i]
        x2 = x

        x0, y0 = newton(y_inicial=y_ini, dominio_x=(x0, x0), n_itens=1)
        x1, y1 = newton(y_inicial=y_ini, dominio_x=(x1, x1), n_itens=1)
        x2, y2 = newton(y_inicial=y_ini, dominio_x=(x2, x2), n_itens=1)

        Dhf = ( y0 - 4*y1 + 3*y2 ) / (2 * vet_h[i])
        vet_dhf[i] = Dhf

    return vet_dhf


def approx_of_df( x0, h ):
    return ( valores_f[ x0 + h ] - valores_f[ x0 - h ] ) / ( 2*h )


def dyF4( y ):
    return 12*( y ** 2 - 2 * y + 2) / (( y - 4 ) ** 3 * y ** 2 * np.sqrt(-(y-4)* y ) )


def max_dF4():
    #first value of dict
    max_ = np.abs( dyF4( list( valores_f.values() )[ 0 ]  ) )
    for i in valores_f:
        if np.abs( dyF4( valores_f[ i ] ) ) > max_:
            max_ = np.abs( dyF4( valores_f[ i ] ) )

    print("Maximo da derivada quarta de f é: %.15e" % np.abs( max_ ) )
    return np.abs( max_ )


#descobrindo n para simpson
def define_n_for_simpson( erro_desej = 10e-05):
    max_f = max_dF4()
    h = ( ( 180 * erro_desej ) / ( (_fim - _inicio) * max_f ) ) ** 1/4
    return np.ceil( ( _fim - _inicio ) / h )



def arc_length(x):
    return np.sqrt( 1 + approx_of_df( x ) ** 2 )


# INTEGRAÇÃO
def G():
    pass


def integral(intervalo=(_inicio, _fim), n=501):
    h = (intervalo[1] - intervalo[0]) / n
    soma_j1 = 0
    soma_j2 = 0
    for i in range(1, n-1):
        if i % 2 == 0:
            soma_j1 += G(X[i])
        else:
            soma_j2 += G(X[i])
    return h / 3 * (G(X[0]) + 2 * soma_j1 + 4 * soma_j2 + G(X[n-1]))


X, Y = newton()
#print("X:")
#print(X)
#print("Y:")
#print(Y)
#print()


#dicionario mapeando domínio e contra domínio
valores_f = {}
for i in range(len(X)):
    valores_f[X[i]] = Y[i]

print("Valor extremo x = %.2f: %.15e" % (_inicio, valores_f[_inicio]))
print("Valor médio x = %.2f: %.15e" % (_meio, valores_f[_meio]))
print("Valor extremo x = %.2f: %.15e" % (_fim, valores_f[_fim]))

#plt.plot(X, Y)
# plt.grid()
# plt.show()


# Derivadas, não é pra fazer daquele jeito que estávamos fazendo
# É pra fazer várias vezes para o ponto e então decidir o tamanho do passo e o tamanho do h
# que nem o professor fez no lab (provavelmente na folha do lab tem explicando isso)

vet_h = np.zeros((10,))
for i in range(len(vet_h)):
    vet_h[i] = np.power(10.0, -(i+2))
    #vet_h = [ 0.1, 0.01, 0.001, 0.0001, ... ]


#derivada de três pontos calculada no ponto central 1
_df1 = df_to_determine_suitable_h_centered( 1 )
print("[~]")
for i, x in enumerate(_df1):
    print("Para h = %.1e => df(1) = %.15e" % (vet_h[i], x))

df1 = df_to_determine_suitable_h_centered( _meio )
df2 = df_to_determine_suitable_h_forward( _inicio )
df3 = df_to_determine_suitable_h_backward( _fim )

# Disso tiramos que o h entre 10^-6 e 10^-8 é OK.
# Podemos verificar os valores nos prints e no plot baixo
# Não sabemos como descobrir O mais preciso. Perguntar ao professor
#for i, x in enumerate(df1):
#    print("Para h = %.1e => df(1) = %.15e" % (vet_h[i], x))
#h ótimo é 10e-06 para fórmula de 3 pontos centrada
#optimal_h = 10e-06

print("Descobrindo h ótimo\n Meio (centrada):\n")
for i, x in enumerate(df1):
    print("Para h = %.1e => df(%d) = %.15e" % (vet_h[i], _meio, x))

print("\nInicio (avançada):")
for i, x in enumerate(df2):
    print("Para h = %.1e => df(%d) = %.15e" % (vet_h[i], _inicio, x))

print("\nFinal: (atrasada)")
for i, x in enumerate(df3):
    print("Para h = %.1e => df(%d) = %.15e" % (vet_h[i], _fim, x))


print("N para simpson: %.15e" % define_n_for_simpson())


# plt.plot(np.log10(vet_h), np.log10(df1))
# plt.grid()
# plt.show()