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

[Leon Mosquera Rodriguez]

 Hola Luis,
de Python no te digo nada, tengo que leer y empezar como autodidacta, factores externos me vienen atrasando en este frente. Conociendo el tema de viscosidad a nivel fisicoquímico, me parece que con un solo punto (50ºC) no podés determinar los valores a 40 ni 100ºC. Resulta que la viscosidad se modela con una ecuación de dos parámetros, y es del tipo exponencial "tipo Arrhenius": v = v' exp(k/T), donde T es temperatura absoluta, y por fiteo o tanteo deberías sacar v' y k para resolver "v". Dos puntos como mínimo, son el input. Como veo logaritmos creo que el modelo está implícito por lo que te dice la Norma. Con tres puntos deberías hacer una regresión linealizando todo: log(v)=v'+k/T. O sea graficando log(v) contra "1/T" tenés una ordenada al origen (v') y la pendiente k. Con Excel sale más rápido...
Si todo esto es superfluo/conocido para vos, desechalo. En cuanto a combinar dos ecuaciones y buscar raíces, ya con dos variables estamos jodidos. Resulta que incluso si las ecuaciones son más o menos "tratables" podrías encontrar infinitas raíces. Te lo simplifico, agarremos un paraboloide para una función de dos variables z=f(x,y), y lo "bajamos" para que corte el plano xy. Las raíces de z=x**2+y**2-c (con c positivo) son un círculo, o sea infinitos puntos. Si no podés reducir el problema a y=f(x), o sea una ecuación (por horrible que sea) con solo una incógnita, no hay tanteo que valga, cuando operes dos ecuaciones, dos variables y un único punto para pivotear. Necesitás añadir una ecuación que opere como "una restricción" (sin que sea necesario caer en los multiplicadores de Lagrange), que a veces es simplemente meter un punto más.
Disculpen la interrupción, saludos!
                                                                       León
    El martes, 13 de marzo de 2018 16:46:49 ART, Luis Andraschnik <luis.andraschnik en gmail.com> escribió:  
 
 
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!!!



 
|  | Libre de virus. www.avast.com  |

_______________________________________________
Lista de Correo de PyAr - Python Argentina - pyar en python.org.ar
Sitio web: http://www.python.org.ar/

Para administrar la lista (o desuscribirse) entrar a http://listas.python.org.ar/listinfo/pyar

La lista de PyAr esta Hosteada en USLA - Usuarios de Software Libre de Argentina - http://www.usla.org.ar  
------------ próxima parte ------------
Se ha borrado un adjunto en formato HTML...
URL: <http://listas.python.org.ar/pipermail/pyar/attachments/20180313/6c025e21/attachment-0001.html>