如何绘制一个非常复杂的函数的导数

如何解决如何绘制一个非常复杂的函数的导数

这个问题与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)

解决方法

选项 #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)

,

使用有限差分法计算 dU/dt 和 dP/dt 的近似值并绘制所需的导数。