0

I'm studying an article where there are two functions:

$e(x) = (e_0/8)\cdot[(2x^3 + x)\sqrt{(1 + x^2)} − \sinh^{−1}(x)]$

$p(x) = (e_0/24)\cdot[(2x^3 - 3x)\sqrt{(1 + x^2)} + 3\sinh^{−1}(x)]$

The article ask me to find numerically $e(p)$ starting from these two expressions. The article suggests me to use a "root-finding routine" but I have no idea how to implement the code. Can someone please help me? I need, possibly, a more general numerical algorithm. (I'm writing in Python). I tried with pynverse library but it is not sufficiently accurate.

mcp
  • 23
  • 3
  • Just to be clear, you are not interested in the composition of e and p, but in the solution of `y==p(x)` used to evaluate `e(x)`? That is, in more compact form, `e(p^{-1}(y))`? – Lutz Lehmann Sep 16 '21 at 09:11
  • @LutzLehmann you're right, e(p) is a compact form for e(p^{-1}(y)). – mcp Sep 16 '21 at 15:23
  • Then `24*y/e_0` for `y=p(x)` is approximately between `2*x^4` (large `x`) and `8/5*x^5` (smaller `x`). One could probably take the corresponding `x` values as initial values for the rapidly converging secant method. One could also use a pre-computed function table for a reverse look-up. All-in-all the resulting function looks like a straight line `e=c*p` apart from a small section close to zero with a vertical tangent at zero. – Lutz Lehmann Sep 16 '21 at 16:28

1 Answers1

0

I solved my problem with the following code (it uses the bisection root finding):

import numpy as np 


def bisection(func,a,b,N):
    if func(a)*func(b) >= 0:
        # No solutions.
        return None
    a_n = a
    b_n = b
    for n in range(1,N+1):
        m_n = (a_n + b_n)/2
        f_m_n = func(m_n)
        if func(a_n)*f_m_n < 0:
            a_n = a_n
            b_n = m_n
        elif func(b_n)*f_m_n < 0:
            a_n = m_n
            b_n = b_n
        elif f_m_n == 0:
            # Found exact solution.
            return m_n
        else:
            print("No solutions with bisection. Return 'None'.")
            return None
    return (a_n + b_n)/2


def eden (pressure):
          
    e = lambda x : (1/8)*((2*x**3 + x)*(1 + x**2)**(1/2) - np.arcsinh(x))
    p = lambda x : (2/24)*((2*x**3 - 3*x)*(1 + x**2)**(1/2) + 3*np.arcsinh(x)) - pressure

    x_p = bisection(p, -1e3, 1e3, 1000)


    return e(x_p)
mcp
  • 23
  • 3
  • Bisection is always the worst solution. Take regula falsi variant illinoise, or Dekker's method if you want to code yourself, `scipy.optimize.find_root` for a more sophisticated and well-tested method. – Lutz Lehmann Sep 16 '21 at 15:37
  • The leading factors seem not to conform to the pattern in the question. – Lutz Lehmann Sep 16 '21 at 16:19