In [2]:
import scipy.optimize as opt
import numpy as np
def half_circle(x):
return (1-x**2)**0.5
下面的程序使用经典的分小矩形计算面积总和的方式，计算出单位半圆的面积：
In [3]:
N = 10000
x = np.linspace(-1, 1, N)
dx = 2.0/N
y = half_circle(x)
面积的两倍
In [4]:
dx * np.sum(y[:-1] + y[1:])
Out[4]:
np.float64(3.1412751679989044)
利用上述方式计算出的圆上一系列点的坐标，还可以用 numpy.trapz 进行数值积分:
In [5]:
import numpy as np
面积的两倍
In [6]:
np.trapz(y, x) * 2
/tmp/ipykernel_177/1721417065.py:1: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`. np.trapz(y, x) * 2
Out[6]:
np.float64(3.1415893269315975)
此函数计算的是以x,y为顶点坐标的折线与X轴所夹的面积。
同样的分割点数， trapz 函数的结果更加接近精确值一些。
如果我们调用 scipy.integrate 库中的 quad 函数的话，将会得到非常精确的结果:
In [7]:
from scipy import integrate
pi_half, err = integrate.quad(half_circle, -1, 1)
pi_half*2
Out[7]:
3.1415926535897967
In [8]:
def half_sphere(x, y):
return (1-x**2-y**2)**0.5
X-Y轴平面与此球体的交线为一个单位圆，因此积分区间为此单位圆，
可以考虑为X轴坐标从-1到1进行积分，
而Y轴从 -half_circle(x) 到 half_circle(x) 进行积分，于是可以调用 dblquad 函数:
In [9]:
integrate.dblquad(half_sphere, -1, 1, lambda x:-half_circle(x),lambda x:half_circle(x))
Out[9]:
(2.0943951023931984, 1.0002354500215915e-09)
通过球体体积公式计算的半球体积：
In [10]:
np.pi*4/3/2
Out[10]:
2.0943951023931953
dblquad 函数的调用方式为：
dblquad(func2d, a, b, gfun, hfun)对于 func2d(x,y) 函数进行二重积分，其中a,b为变量x的积分区间，
而 gfun(x) 到 hfun(x) 为变量y的积分区间。
X轴的积分区间为-1.0到1.0,对于 X=x0 时，通过对Y轴的积分计算出切面的面积，
因此Y轴的积分区间如图中红色点线所示。