How to plot the derivative of a very complicated function

user_7187593 0
Viewed 78

This question is related with the Saha equation. So as you can see, the equation is super long and I defined a lot of variables to simplify this long equation. I want to plot (dU/dt)*(dP/dt)^-1 where U and P are all defined in the code. This equation has a lot of variables which makes it complicated to do this task. I have tried all the methods I could and I still cannot get a good result. Please help me, thank you so much.

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)
derivative python
2 Answers

Option #1:
Calculate the derivative analytically, e.g. with Wolframalpha or by hand and plot it.

Option #2:
Use numerical gradient approximations, for example np.gradient(). Here is an example for your problem:

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)

enter image description here

user_112208840

Use the finite difference method to calculate approximations of dU/dt and dP/dt and plot your required derivative.

user_6469960
Related