L.R.
L.R.

Reputation: 445

Solving a mathematical function in Python

I am trying to solve a logarithmic function using Python. I am searching for an irrational number, so I wrote this bisection algorithm:

def racine(x):
   a=0
   b=x/2
   c=(a+b)/2
   while a!=b and a<c<b:
       if c**2<x:           
           a=c
       else:
           b=c            
       c=(a+b)/2        
   return a, b

which seems to work, at least for finding irrational roots. However then I have a more complicated function:

ln(P)=A+B/T+C*ln(T)

where P, A, B and C are known constants. Isolating T, there is this:

T==e**((ln(P)-A-B/T)/C)

But this still can't be solved because T is on both sides. Can somebody see the way around it? For now I have this code, which clearly doesn't work.

def temperature(P): 
   A=18.19
   B=-23180
   C=-0.8858
   T==e**((log(P)-A-B/T)/C)
   return racine (T)

Thank you!

Upvotes: 1

Views: 450

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

The answer should be to use the bisection method again.

a=small estimate
fa = f(a)
b=large estimate
fb = f(b)
while( b-a > 1e-12 ) {
  c = (a+b)/2
  fc = f(c)
  if( fabs(fc) < 1e-12) return c;
  if( (fc>0) == (fa>0) ) {
    a = c; fa = fc
  } else {
    b = c; f = fc;
}
return (a+b)/2

For more efficient methods look up the regula falsi method in its Illinois variant.

Upvotes: 1

xnx
xnx

Reputation: 25528

If you have NumPy installed, you can find the temperature at a given pressure numerically, for example with scipy.optimize.newton. For example,

import numpy as np
from scipy.optimize import newton
A, B, C = 18.19, -23180, -0.8858
fr = lambda T, lnp: (A + B/T + C*np.log(T)) - lnp
def T(p):
    return newton(fr, 1000, args=(np.log(p),))

In [1]: p1 = 10
In [2]: T1 = T(p1)
In [3]: T1
Out[3]: 2597.8167133280913
In [4]: np.exp(A + B/T1 + C*np.log(T1))    # check
Out[4]: 10.000000000000002

The initial guess value (here 1000) you might have to customize for your usage: I don't know your units.

Upvotes: 0

Related Questions