matplotlib 从Mathcad到Python的翻译

szqfcxe2  于 2022-11-15  发布在  Python
关注(0)|答案(1)|浏览(122)

我无法将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()

ee7vknir

ee7vknir1#

如您所见,MathCad将:
1.创建一个包含S1的符号派生的表达式。
1.找到表达式的根。
在Python中,我们必须使用不同的库来实现相同的结果。在这个特殊的例子中,它稍微复杂一些(需要更多的步骤)。特别是:
1.使用SymPy创建一个符号表达式,然后我们可以计算符号导数。
1.使用SciPy的求根算法,例如rootbisect,...
代码如下:我添加了一些注解以帮助您理解。

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

相关问题