Programa de bisseção

[Lucas Ferreira Loz]

Estou fazendo um código para calcular o ponto de bolha de uma mistura binaria, usando o método de bisseção. Estou com alguns problemas, esse código demora pra pedir as entradas e não retorna nenhuma resposta. Alguém poderia me dizer onde estou errando?

#calculo da temperatura do ponto de bolha
#função objetivo: F=0.55*(exp(14.5463-(2940.46/T-237.22)))+0.45*(exp(14.2724-(2945.47/T-224.0)))-65
import math
a=int(input('a:'))
b=int(input('b:'))
erro=10
while abs(erro)>0.1:
    c=(a+b)/2
    y1=(760 - 0.5*(math.exp(-3848.0840/a + 17.5318) + math.exp(-4328.1150/a + 17.9134)))
    y2=(760 - 0.5*(math.exp(-3848.0840/c + 17.5318) + math.exp(-4328.1150/c + 17.9134)))
    if y2==0:
        print(c)
    elif y1*y2<0:
        b=c
    else:
        a=c
    erro=y2
print('O valor da raiz é: {:.2f}'.format(c))
 

[ArteEN]

ola, tenho duas sugestões

a primeira é você prover uma entrada valida e o resultado esperado, com isso eu posso tentar ajudar mais doque minha proxima sugestão

a segunda sugestão é você colocar um print(y1, y2) para ver o resultada a cada inteção do loop, assim você vai acompanhar passo a passo os resultado matematicos

[Lucas Ferreira Loz]

para entradas use 80 e 100

e o valor esperado é 365

[ArteEN]

F=0.55*(exp(14.5463-(2940.46/T-237.22)))+0.45*(exp(14.2724-(2945.47/T-224.0)))-65

y1=(760 - 0.5*(math.exp(-3848.0840/a + 17.5318)) + math.exp(-4328.1150/a + 17.9134))

bem, ainda não achei nada de errado no codigo, se poder me exclarecer o porque da diferença da função comentada no inicio do script com a função aplicada no codigo

Editado Dezembro 18, 2018 por ArteEN

[Lucas Ferreira Loz]

isso foi um erro meu. Depois eu mudei a função e esqueci de mudar na legenda

[ArteEN]

bem acho que deve ser algo na formula, mas não tenho certesa, vou colocar aqui a saida de 20 inteirações, é um print da numeração da inteção, seguido do erro, y1 , e y2

 0 | 759.99999999999442934495 | 760.00000000000000000000 | 759.99999999999442934495
 1 | 759.99999999994690824678 | 759.99999999999442934495 | 759.99999999994690824678
 2 | 759.99999999984970600053 | 759.99999999994690824678 | 759.99999999984970600053
 3 | 759.99999999975204900693 | 759.99999999984970600053 | 759.99999999975204900693
 4 | 759.99999999968315478327 | 759.99999999975204900693 | 759.99999999968315478327
 5 | 759.99999999964222752169 | 759.99999999968315478327 | 759.99999999964222752169
 6 | 759.99999999961983121466 | 759.99999999964222752169 | 759.99999999961983121466
 7 | 759.99999999960823515721 | 759.99999999961983121466 | 759.99999999960823515721
 8 | 759.99999999960232344165 | 759.99999999960823515721 | 759.99999999960232344165
 9 | 759.99999999959925389703 | 759.99999999960232344165 | 759.99999999959925389703
10 | 759.99999999959777596814 | 759.99999999959925389703 | 759.99999999959777596814
11 | 759.99999999959698016028 | 759.99999999959777596814 | 759.99999999959698016028
12 | 759.99999999959663909976 | 759.99999999959698016028 | 759.99999999959663909976
13 | 759.99999999959641172609 | 759.99999999959663909976 | 759.99999999959641172609
14 | 759.99999999959629803925 | 759.99999999959641172609 | 759.99999999959629803925
15 | 759.99999999959629803925 | 759.99999999959629803925 | 759.99999999959629803925
16 | 759.99999999959629803925 | 759.99999999959629803925 | 759.99999999959629803925
17 | 759.99999999959629803925 | 759.99999999959629803925 | 759.99999999959629803925
18 | 759.99999999959629803925 | 759.99999999959629803925 | 759.99999999959629803925
19 | 759.99999999959629803925 | 759.99999999959629803925 | 759.99999999959629803925
20 | 759.99999999959629803925 | 759.99999999959629803925 | 759.99999999959629803925

em seguida o codigo com as alterações para exibir desta forma que printei em cima,

 

#calculo da temperatura do ponto de bolha
#função objetivo: F=0.55*(exp(14.5463-(2940.46/T-237.22)))+0.45*(exp(14.2724-(2945.47/T-224.0)))-65
import math
a=80 #int(input('a:'))
b=100 #int(input('b:'))
erro=10
cont = 0
while abs(erro)>0.1:
    c=(a+b)/2
    y1=(760 - 0.5*(math.exp(-3848.0840/a + 17.5318) + math.exp(-4328.1150/a + 17.9134)))
    y2=(760 - 0.5*(math.exp(-3848.0840/c + 17.5318) + math.exp(-4328.1150/c + 17.9134)))
    if y2==0:
        print(c)
    elif y1*y2<0:
        b=c
    else:
        a=c
    erro=y2


    print('{:2d} | {:.20f} | {:.20f} | {:.20f}'.format(cont, abs(erro), y1, y2))

    if cont == 20:
        break

    cont+=1
print('O valor da raiz é: {:.2f}'.format(c))

 

[ArteEN]

bem, eu estive procurando material para tentar ajudar, achei alguns exemplos ate mesmo um feito em c#

https://www.devmedia.com.br/resolucao-de-equacoes-metodo-da-bisseccao/24086

tem o algorítimo e tudo

 

estive analizando a equação e quebrei ela em tres lambdas, talvez fique mais facil compreender

 

f1 = lambda x :(
    exp( -3848.0840 / x + 17.5318 )
    )

f2 = lambda x :(
    exp( -4328.1150 / x + 17.9134 )
    )

f3 = lambda x :(
    760 - 0.5 * ( f1(x) + f2(x) )
    )

 

[ArteEN]

ola novamente, acho que estou chegando em algum lugar aqui, vou citar minha observação

usando a função lambda que apresentei a cima

f1 = lambda x : (
    exp(-3848.0840/x + 17.5318)
    )

ela retorna um numero minusculo, você pode conferir com este codigo

c = f1(a)
d = 0.0


for i in range(10):
    print( c > d )

    d += 0.0000000000001

 

você pode ver que estou incrementando a comparação esta em menos do que 1 milionésimo, isso deve explicar do porque esta demorando tanto, o ponto é a formula que você esta usando, ela esta numa precisão muito grande para se ter uma computação rapida

 

só para você ter uma ideia a saida da formula descartando o 760 é de 2.6575817451649962e-14

repare no final do numero e-14 isso é deslocamento de virgula, para ver o numero sem o deslocamento eu usei um print formatado, eis o resultado 0.00000000000002657581745164996181812582192510641086

com essa observação eu apliquei o print formatado por dois ciclos do loop, e você pode ver o resultado a seguir

 

1 = 759.99999999999442934495164081454277038574218750000000
2 = 759.99999999994690824678400531411170959472656250000000

 

de um ciclo para o outro a mudança é minima, e é por isso a demora