[Python] Python Científico: Código para calcular a temperatura do ponto de Bolha usando 3...

estou tentando criar um código que calcula a temperatura do ponto de bolha para a gasolina, os dados foram:

A gasolina pode ser aproximada por uma mistura de 4 alcanos: C6: 15% em mol C7: 25% em mol C8: 45% em molC9: 15 % em mol

Equação de Antoine para as pressões de saturação: P_i, sat = exp [A_i - B_i / (T + C_i)], em que P_i,sat [=] mmHg e T [=] K C6: A = 15.8366; B = 2697,50; C= -48,78 C7: A = 15.8737; B = 2911,32; C = -56,51 C8: A = 15,9426; B = 3120,29; C = -63,63 C9: A = 15,9671; B = 3291,45; C = -71,33

e eu tenho que fazer o código usando 4 métodos diferentes para 3 funções diferentes que estão no código, o problema acaba sendo na última função, vocês poderiam me ajudar a deixar o código mais otimizado para que eu consiga calcular a última função?

import math as mt
import numpy as np
import matplotlib.pyplot as plt

'Author: @nogfe4'

# 1: a) Newton-Rapshon

# Pressão:
P = 760 # [mmhg] ou 1 [atm]

# Dados da Gasolina (Mistura de Hidrocarbonetos):
xC6 = 0.15
AC6 = 15.8366
BC6 = 2697.50
CC6 = -48.78
xC7 = 0.25
AC7 = 15.8737
BC7 = 2911.32
CC7 = -56.51
xC8 = 0.45
AC8 = 15.9426
BC8 = 3120.29
CC8 = -63.63
xC9 = 0.15
AC9 = 15.9671
BC9 = 3291.45
CC9 = -71.33

def f(T):
f = ((np.exp(AC6 - BC6/(T+CC6))*xC6) + ((np.exp(AC7-BC7/(T+CC7)))*xC7)
+ (np.exp(AC8-BC8/(T+CC8))*xC8) + (np.exp(AC9-BC9/(T+CC9))*xC9)) - P
return f

# Calculando o chute a partir das temperaturas de saturação:
T6_sat = BC6/(AC6 - np.log(P)) - CC6
T7_sat = BC7/(AC7 - np.log(P)) - CC7
T8_sat = BC8/(AC8 - np.log(P)) - CC8
T9_sat = BC9/(AC9 - np.log(P)) - CC9

T0 = (T6_sat*xC6 + T7_sat*xC7 + T8_sat*xC8 + T9_sat*xC9)

# Dados da iteração
it = 0
itmax = 100
tol = 1e-6

# Derivada numérica da Função:
def numderivative(f, T, h=1e-6):
nd = (f(T + h) - f(T - h)) / (2 * h)
return nd

# 1: a) Newton-Rapshon
# Cálculo da iteração:
while it < itmax:
f0 = f(T0)
df0 = numderivative(f, T0)
T_new = T0 - f0/df0

if abs(T_new - T0) < tol:
break
T0 = T_new
it = it + 1

# Resultados:
print("************************** EXERCÍCIO 1 ***********************************************************")
print(" Exercício 1 - Usando o método de Newton-Raphson")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(T_new, 2)} K")
print("---------------------------------------------------------------------------------------------------")



# 1: b) Secante
# Os dados e as funções são as mesmas, oq muda é o algoritmo do método e a inserção de uma nova temperatura T1

T1 = T0*0.98
tol2 = 1e-6
# Cálculo da iteração:
while it < itmax:
f0 = f(T0)
f1 = f(T1)
T_new2 = T0 - ((f0*(T1-T0))/(f1-f0))

if abs(T_new2 - T0) < tol2:
break
it +=1

print(" Exercício 1 - Usando o método da Secante")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(T_new2, 2)} K")
print("---------------------------------------------------------------------------------------------------")

# 1: c) Regula Falsi
# Os dados e as funções são as mesmas, oq muda é o algoritmo do método e a inserção de novos dados


# Cálculo da iteração:
# Supondo que f(T), Ti, Tf, it, itmax e tol já foram definidos

Ti = 298
Tf = 400

while it < itmax:
T_new3 = Tf - ((f(Tf)*(Ti-Tf))/(f(Ti)-f(Tf)))

if abs(f(T_new3)) < tol:
break
if f(Ti)*f(T_new3) < 0:
Tf = T_new3
else:
Ti = T_new3
it = it + 1

if it == itmax:
print("O método não convergiu no número máximo de iterações.")



print(" Exercício 1 - Usando o método da Falsa Posição (Regula Falsi)")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(T_new3, 2)} K")
print("---------------------------------------------------------------------------------------------------")


# 1: d) Bisseção
# Os dados e as funções são as mesmas, oq muda é o algoritmo do método e a inserção de novos dados

Ti = 300
fi = f(Ti)
Tf = 391
ff = f(Tf)

if f(Ti) * f(Tf) >= 0:
print("Intervalo inicial para Bisseção é inválido")

else:
Tc = (Ti + Tf)/2

while it < itmax:
Tc = (Ti + Tf)/2
fc = f(Tc)
if abs(fc) < tol:
break

if fi*fc < 0:
Tf = Tc
else:
Ti = Tc
fi = fc
it = it + 1


print(" Exercício 1 - Usando o método da Bisseção")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(Tc, 2)} K")
print("---------------------------------------------------------------------------------------------------")


# 2: a) Newton Raphson
# Definindo a nova função f2 que será usada no item 2

def f2(T):
f2 = np.log(((np.exp(AC6 - BC6/(T+CC6))*xC6) + ((np.exp(AC7-BC7/(T+CC7)))*xC7)
+ (np.exp(AC8-BC8/(T+CC8))*xC8) + (np.exp(AC9-BC9/(T+CC9))*xC9))) - np.log(P)
return f2


# Derivada numérica da Função:
def numderivative(f2, T, h=1e-6):
nd = (f2(T + h) - f2(T - h)) / (2 * h)
return nd

# Cálculo da iteração:
while it < itmax:
f0 = f2(T0)
df0 = numderivative(f2, T0)
T_new4 = T0 - f0/df0

if abs(T_new4 - T0) < tol:
break
T0 = T_new4
it = it + 1

# Resultados:
print("************************** EXERCÍCIO 2 ***********************************************************")
print(" Exercício 2 - Usando o método de Newton-Raphson")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(T_new4, 2)} K")
print("---------------------------------------------------------------------------------------------------")

# 2: b) Secante

T1 = T0*0.98
tol2 = 1e-6

# Cálculo da iteração:
while it < itmax:
f0 = f2(T0)
f1 = f2(T1)
T_new5 = T0 - ((f0*(T1-T0))/(f1-f0))

if abs(T_new5 - T0) < tol2:
break
it +=1

print(" Exercício 2 - Usando o método da Secante")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(T_new5, 2)} K")
print("---------------------------------------------------------------------------------------------------")

# 2: c) Regula Falsi
# Os dados e as funções são as mesmas, oq muda é o algoritmo do método e a inserção de novos dados


# Cálculo da iteração:
# Supondo que f(T), Ti, Tf, it, itmax e tol já foram definidos

Ti = 298
Tf = 400

while it < itmax:
T_new6 = Tf - ((f2(Tf)*(Ti-Tf))/(f2(Ti)-f2(Tf)))

if abs(f2(T_new6)) < tol:
break
if f2(Ti)*f2(T_new6) < 0:
Tf = T_new6
else:
Ti = T_new6
it = it + 1

if it == itmax:
print("O método não convergiu no número máximo de iterações.")



print(" Exercício 2 - Usando o método da Falsa Posição (Regula Falsi)")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(T_new6, 2)} K")
print("---------------------------------------------------------------------------------------------------")


# 2: d) Bisseção
# Os dados e as funções são as mesmas, oq muda é o algoritmo do método e a inserção de novos dados

Ti = 300
fi = f2(Ti)
Tf = 391
ff = f2(Tf)

if f2(Ti) * f2(Tf) >= 0:
print("Intervalo inicial para Bisseção é inválido")

else:
Tc2 = (Ti + Tf)/2

while it < itmax:
Tc2 = (Ti + Tf)/2
fc = f2(Tc2)
if abs(fc) < tol:
break

if fi*fc < 0:
Tf = Tc2
else:
Ti = Tc2
fi = fc
it = it + 1


print(" Exercício 2 - Usando o método da Bisseção")
print("---------------------------------------------------------------------------------------------------")
print(f"Temperatura do ponto de Bolha para a Gasolina à pressão atmosférica (760 mmHg) é: {round(Tc, 2)} K")
print("---------------------------------------------------------------------------------------------------")



# 3: a) Newton-Rapshon
# Definindo a nova função f3 que será usada no item 2

j = 1/T0

def f3(j):
f = math.log(((AC6 - BC6)/(T+CC6)*xC6) + (AC7-BC7/(T+CC7)*xC7)
+ (AC8-BC8)/(T+CC8)*xC8) + (AC9-BC9)/(T+CC9))*xC9)) - math.log(P)
return f

def f3_derivada(j, h=1e-6):
return (f3(j + h) - f3(j - h)) / (2*h)


Onde posso melhorar?

Continue reading...