python - 如何绘制非常复杂函数的导数

Question

此问题与 Saha equation 相关。正如您所看到的，方程非常长，我定义了很多变量来简化这个长方程。我想绘制 (dU/dt)*(dP/dt)^-1 ，其中 U 和 P 均在代码中定义。该方程有很多变量，这使得完成这项任务变得复杂。我已经尝试了所有能用的方法，但仍然得不到好的结果。请帮助我，非常感谢。

import numpy as np
import math
import matplotlib.pyplot as plt
import sympy as sym  

# These are the constant variables I defined 
kb = 1.381e-23
h = 6.626e-34
me = 9.109e-31
mh = 1.67e-27
pi = math.pi
rho1 = 10e-7
rho2 = 10e-5
rho3 = 10e-3


#the following variables will be used in the definition of P and U
t = np.linspace(5000, 25000, 1000)
k1 =(1/rho1)*(mh*((2*3.1415926*me*kb)**(3/2)))/(h**3)
k2 = 2.18e-18/kb 
k3 =(1/rho2)*(mh*((2*3.1415926*me*kb)**(3/2)))/(h**3)
k4 =(1/rho3)*(mh*((2*3.1415926*me*kb)**(3/2)))/(h**3)

x1 = ((-k1* (t**(3/2))*np.exp(-k2/t)+np.sqrt((k1**2)*(t**3)*np.exp(-k2*2/t)+4*k1*(t**(3/2))*np.exp(-k2/t)))/2)
x2 = ((-k3* (t**(3/2))*np.exp(-k2/t)+np.sqrt((k3**2)*(t**3)*np.exp(-k2*2/t)+4*k3*(t**(3/2))*np.exp(-k2/t)))/2)
x3 = ((-k4* (t**(3/2))*np.exp(-k2/t)+np.sqrt((k4**2)*(t**3)*np.exp(-k2*2/t)+4*k4*(t**(3/2))*np.exp(-k2/t)))/2)

# I defined P and U here 
P = (rho3/mh)*(1+x3)*(kb*t)
U = 1.5*P + rho3*x3*2.18e-18/(mh)

Answer 1

选项#1:
分析计算导数，例如与 Wolframalpha或者手动绘制它。

选项#2:
使用数值梯度近似，例如 np.gradient() 。以下是您的问题的示例:

def saha(t):
    kb = 1.381e-23
    h = 6.626e-34
    me = 9.109e-31
    mh = 1.67e-27
    pi = math.pi
    rho1 = 10e-7
    rho2 = 10e-5
    rho3 = 10e-3
    k1 =(1/rho1)*(mh*((2*3.1415926*me*kb)**(3/2)))/(h**3)
    k2 = 2.18e-18/kb 
    k3 =(1/rho2)*(mh*((2*3.1415926*me*kb)**(3/2)))/(h**3)
    k4 =(1/rho3)*(mh*((2*3.1415926*me*kb)**(3/2)))/(h**3)

    x1 = ((-k1* (t**(3/2))*np.exp(-k2/t)+np.sqrt((k1**2)*(t**3)*np.exp(-k2*2/t)+4*k1*(t**(3/2))*np.exp(-k2/t)))/2)
    x2 = ((-k3* (t**(3/2))*np.exp(-k2/t)+np.sqrt((k3**2)*(t**3)*np.exp(-k2*2/t)+4*k3*(t**(3/2))*np.exp(-k2/t)))/2)
    x3 = ((-k4* (t**(3/2))*np.exp(-k2/t)+np.sqrt((k4**2)*(t**3)*np.exp(-k2*2/t)+4*k4*(t**(3/2))*np.exp(-k2/t)))/2)
    P = (rho3/mh)*(1+x3)*(kb*t)
    U = 1.5*p + rho3*x3*2.18e-18/(mh)
    return P, U

# compute U and P
t = np.linspace(5000, 25000, 1000)
p, u = saha(t)

# compute gradients of U and P w.r.t. T
u_grad = np.gradient(u)
p_grad = np.gradient(p)

# plot
plt.plot(t,u_grad * 1/(p_grad))
plt.ylabel("(dU/dt)*(dP/dt)^-1")
plt.xlabel("T")
plt.savefig("saha.png",dpi=100)