import numpy as np from scipy import integrate def half_circle(x): return (1 - x ** 2) ** 0.5 N = 10000 x = np.linspace(-1, 1, N) dx = 2. / N; y = half_circle(x) area =sum( dx * y)#利用矩形面积法 print np.trapz(y, x) * 2#求数值积分 pi_half, err = integrate.quad(half_circle, -1,1) #求积分 print pi_half * 2 def half_sphere(x ,y): return (1 - x ** 2 - y ** 2) ** 0.5 print integrate.dblquad(half_sphere,-1,1,lambda x: -half_circle(x),lambda x:half_circle(x))#求二重积分 from scipy.integrate import odeint def lorenz(w ,t, p, r, b): x ,y, z = w return np.array([p * (y -x), x * (r-z)-y, x * y - b * z]) t = np.arange(0 , 40, 0.01) track1 = odeint(lorenz, (0.0,1.00,0.0),t, args=(10.0, 28.0,3.0)) track2 = odeint(lorenz, (0.0,1.01,0.0),t, args=(10.0, 28.0,3.0)) from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt fig = plt.figure() #ax = fig.gca(projection = '3d') ax = Axes3D(fig) ax.plot(track1[:,0], track1[:,1], track1[:,2]) ax.plot(track2[:,0], track2[:,1], track2[:,2]) plt.show()