从Mathcad到Python的翻译

Question

我不能把公式从Mathcad翻译成Python。卡在"a“上。以下是我所能做的：

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

如您所见，MathCad将：

1. 创建一个包含S1符号导数的表达式。
2. 找到那个表达式的根。

在Python中，我们必须使用不同的库来实现相同的结果。在这个特殊的例子中，它更复杂一些(需要更多的步骤)。特别是：

1. 使用SymPy创建一个符号表达式。然后我们可以计算符号导数。
2. 使用SciPy的寻根算法，如root或bisect，.

下面是代码:我添加了一些注释来帮助您理解。

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()

这是以下代码的输出：

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)

我们想要找出根的函数

​

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

​

​