[pyar] Resolviendo ecuaciones con el método de bisección : Ayuda matemátcos pythonistas!!

[Luis Andraschnik]

Estoy tratando de resolver estas ecuaciones para reemplazar unos programas
en DOS muy viejos:

#########Cálculo de índice de viscosidad VI (viscosity index, variación de
viscosidad con la temperatura entre 40 y 100°C):

from bisect import bisect

from math import log, log10

def calc_VI(kv100, kv40):
    """"Calcula índice de viscosidad según ASTM D 2270
        kv100 viscosidad a 100°C en mm2/s (cSt),
        kv40 viscosidad a 40°C en mm2/s (cSt)"""

    kv100_range = (2.0, 3.8, 4.4, 5.0, 6.4, 7.0,
                   7.7, 9.0, 12, 15, 18, 22, 28,
                   40, 55, 70)

    LH_data = ((None),
               (1.14673, 1.7576, -0.109, 0.84155, 1.5521, -0.077),
               (3.38095, -15.4952, 33.196, 0.78571, 1.7929, -0.183),
               (2.5000, -7.2143, 13.812, 0.82143, 1.5679, 0.119),
               (0.10100, 16.6350, -45.469, 0.04985, 9.1613, -18.557),
               (3.35714, -23.5643, 78.466, 0.22619, 7.7369, -16.656),
               (0.01191, 21.4750, -72.870, 0.79762, -0.7321, 14.610),
               (0.41858, 16.1558, -56.040, 0.05794, 10.5156, -28.240),
               (0.88779, 7.5527, -16.600, 0.26665, 6.7015, -10.810),
               (0.76720, 10.7972, -38.180, 0.20073, 8.4658, -22.490),
               (0.97305, 5.3135, -2.200, 0.28889, 5.9741, -4.930),
               (0.97256, 5.2500, -0.980, 0.24504, 7.4160, -16.730),
               (0.91413, 7.4759, -21.820, 0.20323, 9.1267, -34.230),
               (0.87031, 9.7157, -50.770, 0.18411, 10.1015, -46.750),
               (0.84703, 12.6752, -133.310, 0.17029, 11.4866, -80.620),
               (0.85921, 11.1009, -83.19, 0.17130, 11.3680, -76.940),
               (0.83531, 14.6731, -216.246, 0.16841, 11.8493,-96.947)
              )

    U, Y = kv40, kv100

    if Y > 70:
        L = 0.8353*Y*Y + 14.67*Y - 216
        H = 0.1684*Y*Y + 11.85*Y - 97
    else:
        i = bisect(kv100_range, Y)
        a, b, c, d, e, f = LH_data[i]
        L = a*Y*Y + b*Y + c
        H = d*Y*Y + e*Y + f

    VI = (L - U)*100/(L - H)

    if VI > 100:
        N = (log(H) - log(U))/ log(Y)
        VI = (10**N - 1)/0.00715 + 100

    return VI


#########Cálculo de viscosidad vs temperatura:

def viscosidad_a_t3(v1, v2, t1 , t2 , t3):  ########### Ecuación 1
    """v1,v2 viscosidades en mm2/s a t1, t2, temperaturas en °C
       v3 viscosidad requerida para la temperatura t3
       según ASTM D 341"""

    t1+=273.15
    t2+=273.15
    t3+=273.15
    LT1 = log10(t1)
    LT2 = log10(t2)
    LT3 = log10(t3)
    LLv1 = log10(log10(v1 + 0.7))
    LLv2 = log10(log10(v2 + 0.7))
    m = (LLv1 - LLv2) / (LT1 - LT2)
    LLv3 = LLv1 - m*(LT1 - LT3)
    v3 = 10**10**LLv3 - 0.7
    return v3


#####En el primer caso para calcular VI requiere kv40 y kv100 (datos
experimentales). el cálculo es explícito.

##### Para calcular kv40 ó  kv100 conociendo el VI utilizo unas ecuaciones
auxiliares para poder aplicar el método de bisección


def func_kv100(kv100, kv40, iv=70):
    return iv - calc_VI(kv100, kv40)

def func_kv40(kv40, kv100, iv=70):
    return iv - calc_VI(kv100, kv40)

y aplico el algoritmo de bisección:

def biseccion(f, a, b, tol=1.0e-6):
    """Método de bisección
    Halla una raíz de la función f en el intervalo [a, b] mediante el método
    de bisección.....

    https://www.pybonacci.org/2012/04/18/
    ecuaciones-no-lineales-metodo-de-biseccion-y-metodo-de-newton-en-python/
    """
    if a > b:
        raise ValueError("Intervalo mal definido")
    if f(a) * f(b) >= 0.0:
        raise ValueError("La función debe cambiar de signo en el intervalo")
    if tol <= 0:
        raise ValueError("La cota de error debe ser un número positivo")
    x = (a + b) / 2.0
    while True:
        if b - a < tol:
            return x
        # Utilizamos la función signo para evitar errores de precisión
        elif np.sign(f(a)) * np.sign(f(x)) > 0:
            a = x
        else:
            b = x
        x = (a + b) / 2.0

from functools import partial

    #Calcular kv100 con datos kv40 e iv
    f = partial(func_kv100,kv40=kv40, iv=iv)
    kv100 = biseccion(f, a=2, b=100, tol=0.001)  ############### Ecuación 2

    #Calcular kv40 con datos kv100 e iv
    f = partial(func_kv40,kv100=kv100, iv=iv)
    kv40 = biseccion(f, a=2, b=100, tol=0.001)

 ... y si llegaron hasta acá, mi gran problema es el siguiente. Si conozco
la viscosidad a  por ejemplo 50°C y el IV, cómo resuelvo la viscosidad a 40
y 100°C.

Si lo quiero hacer tanteando empiezo proponiendo una viscosidad a 40°C y
calculo la viscosidad a 100 (ecuación 2)

Con la viscosidad a 40°C y 100°C  calculo la viscosidad a 50°C (ecuación 1)

Comparo el valor calculado de viscosidad a 50°C con el dato experimental y
según sea mayor o menor "retoco un poco el valor de la viscosidad a  40°C y
vuelvo hasta alcanzar la presición requerida.

El problema es que no se me ocurre como combinar las 2 ecuaciones para
aplicar el método de bisección ....

Ayuda matemáticos  pythonistas!!!


<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail>
Libre
de virus. www.avast.com
<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail>
<#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
------------ próxima parte ------------
Se ha borrado un adjunto en formato HTML...
URL: <http://listas.python.org.ar/pipermail/pyar/attachments/20180313/ffc4ee66/attachment.html>