Translation from Mathcad to Python

Question

I can't translate a formula from Mathcad to Python. Stuck on "a". Here's what I was able to do:

from matplotlib import pyplot as plt
import numpy as np
k1 = 1
b = 1.51
D = (1/b) * (np.sqrt(k1/np.pi))
x0 = 10 * b

myArray = np.arange(0, 24, 0.1)

for t in myArray:
    S1_t = (k1) / (1 + np.e ** (-(D * myArray - 5)))
    S1_deistv = S1_t.real
plt.plot(myArray, S1_deistv, color="black")
plt.show()

Answer 1

As you can see, MathCad is going to:

1. create an expression containing the symbolic derivative of S1.
2. finding the root of that expression.

In Python, we have to use different libraries to achieve the same result. In this particular case it is a little bit more convoluted (requires more steps). In particular:

1. Use SymPy to create a symbolic expression. We can then compute the symbolic derivative.
2. Use SciPy's root finding algorithms, such as root or bisect, ...

Here is the code: I've added some comments to help you understand.

from matplotlib import pyplot as plt
import numpy as np
from scipy.optimize import root
import sympy as sp

k1 = 1
b = 1.51
D = (1/b) * np.sqrt(k1 / np.pi)

# create a symbol
t = sp.symbols("t")
# create a symbolic expression. Note that we are using the
# exponential function of SymPy (because it is symbolic)
S1 = k1 / (1 + sp.exp(-(D * t - 5)))
print("S1 = ", S1)
# write the expression
# with `diff` we are computing the derivative with respect to `t`
expr = S1 - t * sp.diff(S1, t)
print("expr = ", expr)
# convert the expression to a numerical function so that it
# can be evaluated by Numpy/Scipy
f = sp.lambdify([t], expr)
# plot the symbolic expression to help us decide a good initial
# guess for the root finding algorithm
sp.plot(expr, (t, 0, 24))
# in the interval t in [0, 24] we can see two roots, one at
# about 2 and the other at about 18.
# Let's find the second root.
result = root(f, 18)
print("result", result)
a = result.x[0]
print("a = ", a)

# remember that S1 is a symbolic expression: we can substitute
# t (the symbol) with a (the number)
b = float(S1.subs(t, a))
k = b / a
print("k = ", k)

t_array = np.arange(0, 24, 0.1)
plt.figure()
S1_t = (k1) / (1 + np.e ** (-(D * t_array - 5)))
S1_deistv = S1_t.real
plt.plot(t_array, S1_deistv, color="black")
plt.plot(t_array, k * t_array)
plt.show()

This is the output of the following code:

S1 =  1/(exp(5 - 0.373635485793216*t) + 1)
expr =  -0.373635485793216*t*exp(5 - 0.373635485793216*t)/(exp(5 - 0.373635485793216*t) + 1)**2 + 1/(exp(5 - 0.373635485793216*t) + 1)

Function of which we want to find the roots

result     fjac: array([[-1.]])
     fun: array([-6.66133815e-16])
 message: 'The solution converged.'
    nfev: 6
     qtf: array([3.5682568e-13])
       r: array([-0.22395716])
  status: 1
 success: True
       x: array([18.06314347])
a =  18.063143471730815
k =  0.04715849105203411