我想使用和绘制某个函数的数值积分函数。我怎样才能做到这一点？scipymatplotlib
我尝试了以下但它没有工作（运行ipython %pylab）：
import numpy as np
from scipy import integrate
def f(x):
return x*np.sin(1/x)
X = np.arange(-0.5,0.5,0.001)
#plot(X,f(X))
def F(x):
return integrate.quad(f,0,x)
plot(X,F(X))
我还尝试对函数进行矢量化但F没有成功。
首先，您的函数中是奇异的。您可能想要添加这样的子句：if
def f(x):
if abs(x) < 1e-10:
res = x
else:
res = x*sin(1/x)
但这确实会损害速度（屏蔽数组会更好）。
您的代码不起作用的原因是因为
scipy.quad只接受一个值作为左右边界
scipy.quad返回一个(val,err)包含积分值和误差估计值的元组，因此您需要将其解包。
此代码有效：
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
%matplotlib inline
def f(x):
if (np.abs(x)<1e-10):
res = x
else:
res = x*np.sin(1.0/x)
return res
X = np.arange(-0.5,0.5,0.001)
#plot(X,f(X))
def F(x):
res = np.zeros_like(x)
for i,val in enumerate(x):
y,err = integrate.quad(f,0,val)
res[i]=y
return res
plt.plot(X,F(X))
将最后一行替换为
plot(X, [F(x)[0] for x in X])
那应该这样做。
编辑：您可以将函数 F 定义为
def F(x):
try:
return [integrate.quad(f, 0, y) for y in x]
except TypeError:
return integrate.quad(f, 0, x)
from scipy.integrate import cumtrapz
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return x*np.sin(1/x) if abs(x) > 1e-10 else x
f = np.vectorize(f)
X = np.arange(-0.5, .5, 0.001)
fv = f(X)
plt.plot(fv)
F = cumtrapz(fv, x=X, initial=0)
plt.plot(F);
如果您坚持使用quad，则更有效的实现将使用 quad 计算细分段上的积分以获得最佳精度，然后计算每个样本点的反导数值的累积和。
def f(x):
return x*np.sin(1/x)
X = np.arange(-0.5,0.5,0.001)
DF = [ integrate.quad(f,a,b)[0] for a,b in zip(X[:-1],X[1:]) ]
F = np.cumsum(DF)
plot(X[1:],F)
