Search code examples
pythonwolfram-mathematicascientific-computingquad

python plot integration result , how to do like what Mathematica does


I want to rewrite this Mathematica code into Python, but I was puzzled then, please help me! Thanks a lot. Perhaps the result of an integration is different with array, however, I don' t know!

IMAGE: mathematica code which I want to rewrite in python

I was troubled with python code then....

# -*- 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()

Solution

  • I think that's it ! Hallelujah !!

    #-----------  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--' )
    
    # ------------------------------------------