python 绘制积分结果,如何像 Mathematica 那样做
python plot integration result , how to do like what Mathematica does
我想把这段Mathematica代码重写成Python,但是当时我很疑惑,请帮帮我!非常感谢。也许积分的结果与阵列不同,但是,我不知道!
IMAGE: mathematica code which I want to rewrite in python
当时我被 python 代码困扰....
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
import scipy as sp
from pylab import *
#from scipy import *
from scipy.integrate import quad, dblquad, tplquad
plt.figure('God Bless : fig1')
plt.title(r'fig1_')
plt.xlabel(r'z')
plt.ylabel('$L_0$(erg $s^{-1}$)')
#plt.axis([0,10,-2,2])
Limit = 1.e48
c = 2.997e10
Om = 0.27
H0 = 70e5/(1.e6*3.86e18)
z = np.arange(0,10, 0.01)
def dLz(z):
return 1./(1.-Om + Om*(1.+z)**3.)**(1./2)
val, abserr = quad(dLz, 0, 10)
print ("integral value =", val, ", absolute error =", abserr)
dL = val
Fmin = 2.0e-8
Llimit = 4.0*np.pi*dL**2*Fmin
plt.plot(z, Limit)
plt.savefig('fig_1.eps', dpi=300)
plt.show()
我想就是这样!哈利路亚!!
#----------- DETERMIN zimax -------------
Limit = 1.e48
c = 2.997e10
Om = 0.27
H0 = 70e5/(1.e6*3.86e18)
Fmin = 2.0e-8
z = np.arange(0,10, 0.1)
def dLz(z):
return (c/H0)*(1.+z)/(1.-Om + Om*(1.+z)**3.)**(1./2)
x_lower = 0
vals = []
Llimits = []
for x_upper in z :
val, abserr = quad(dLz, x_lower, x_upper)
vals.append(val)
print ("integral value =", val, ", absolute error =", abserr)
dL = val
Llimit = 4.0*np.pi*dL**2*Fmin
Llimits.append(Llimit) # add to array
plt.plot(z, Llimits, 'b--' )
# ------------------------------------------
我想把这段Mathematica代码重写成Python,但是当时我很疑惑,请帮帮我!非常感谢。也许积分的结果与阵列不同,但是,我不知道!
IMAGE: mathematica code which I want to rewrite in python
当时我被 python 代码困扰....
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
import scipy as sp
from pylab import *
#from scipy import *
from scipy.integrate import quad, dblquad, tplquad
plt.figure('God Bless : fig1')
plt.title(r'fig1_')
plt.xlabel(r'z')
plt.ylabel('$L_0$(erg $s^{-1}$)')
#plt.axis([0,10,-2,2])
Limit = 1.e48
c = 2.997e10
Om = 0.27
H0 = 70e5/(1.e6*3.86e18)
z = np.arange(0,10, 0.01)
def dLz(z):
return 1./(1.-Om + Om*(1.+z)**3.)**(1./2)
val, abserr = quad(dLz, 0, 10)
print ("integral value =", val, ", absolute error =", abserr)
dL = val
Fmin = 2.0e-8
Llimit = 4.0*np.pi*dL**2*Fmin
plt.plot(z, Limit)
plt.savefig('fig_1.eps', dpi=300)
plt.show()
我想就是这样!哈利路亚!!
#----------- DETERMIN zimax -------------
Limit = 1.e48
c = 2.997e10
Om = 0.27
H0 = 70e5/(1.e6*3.86e18)
Fmin = 2.0e-8
z = np.arange(0,10, 0.1)
def dLz(z):
return (c/H0)*(1.+z)/(1.-Om + Om*(1.+z)**3.)**(1./2)
x_lower = 0
vals = []
Llimits = []
for x_upper in z :
val, abserr = quad(dLz, x_lower, x_upper)
vals.append(val)
print ("integral value =", val, ", absolute error =", abserr)
dL = val
Llimit = 4.0*np.pi*dL**2*Fmin
Llimits.append(Llimit) # add to array
plt.plot(z, Llimits, 'b--' )
# ------------------------------------------