在 matplotlib 中绘制积分奇异函数时出现的问题

Question

我想在 matplotlib 中绘制（奇异）函数的积分，但我的代码不起作用。从数学上讲，我想要这个：

代码：

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt


def g(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,g(X)) ## Doesn't work

def f(x):
    res = np.zeros_like(x)
    for i,val in enumerate(x):
        y,err = integrate.quad(g,0,val)
        res[i]=y
    return res

plot(X,f(X)) ## Works

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)) ## Doesn't work

（代码是https://scicomp.stackexchange.com/a/21871/9417的改编版本）

所以我无法绘制原始函数 g，因为它表示：

      5 
      6 def g(x):
----> 7     if (np.abs(x)<1e-10):
      8         res = x
      9     else:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

而且我也不能像它所说的那样绘制积分函数：

     17 def f(x):
     18     res = np.zeros_like(x)
---> 19     for i,val in enumerate(x):
     20         y,err = integrate.quad(g,0,val)
     21         res[i]=y

TypeError: 'float' object is not iterable

但是绘制第一个积分 f 是可行的。我怎样才能解决这个问题？在没有运行的情况下，是否有更好的方法将 python 绘制成此类问题？

Answer 1

您的代码中存在一些错误：

当你运行g(X)时，参数是一个数组，而在函数内部你把X当作一个值。 g(x)并不是应用于X的每一个元素，而是应用于整个数组X。这就解释了为什么 The truth value of an array... 错误，因为 x 实际上是整个数组 X。您可以通过两种不同的方式解决它：
1. 将函数 g 映射到 X
  的每个值
```
>>> y = map(X, g)  # y = np.asarray(map(X, g)) if you want a numpy array back
>>> plot(X, y)
```
y = map(X, g) 也可以使用列表推导式写成 y = [g(x) for x in X].
1. 向量化函数g(x)以应用于整个数组
```
def g(X):
    r = X.copy() # create a copy of x
    mask = np.abs(r) < 1e-10   # create a mask with invalid values
    r[mask] *= np.sin(1.0 / r[mask])  # replace those values
    return r
```
但是，如果您选择选项 2，函数 f(x) 将失败，因为它使用单个值调用 g(x)，每个 X 的值调用一次。
您的函数 f(x) 改为使用数组，因为您遍历 X 的每个元素并将 g(x) 应用于每个元素。如果您仍想将 g(x) 应用于 X 的每个元素，请使用 option 1.
您的函数 F(x) 也适用于整个数组，因为您循环遍历其元素。但是，对于您调用 f(x) 的每个元素，它只允许 lists/arrays 作为输入（并且您给出的是一个数字）。我认为 F(x) 是多余的，因为它与 f(x).

根据您的定义，您可以重写方程式如下：

def g(x):
    return x if np.abs(x) < 1e-10 else x * np.sin(1.0/x)

def f(x):
    return integrate.quad(g, 0, x)[0]

def F(x):
    return integrate.quad(f, 0, x)[0]

然后对每个x映射函数得到结果

X = np.arange(0,0.5,0.01)
yg = map(g, X)
yf = map(f, X)
yF = map(F, X)

绘图结果：

import matplotlib.pyplot as plt
X = np.arange(-0.5, 0.5, 0.01)
yf = map(f, X)
plt.plot(X, yf)
plt.show()

一个更快的版本，允许限制积分内的点数：

def g(x):
    return x if np.abs(x) < 1e-10 else x * np.sin(1.0/x)

def f(x, limit=50):
    return integrate.quad(g, 0, x, limit=limit)[0]

def F(x, limit=50):
    return integrate.quad(f, 0, x, args=limit, limit=limit)[0]

然后运行:

X = np.arange(-0.5, 0.5, 0.01)
FY = [F(x, limit=10) for x in X]
plt.plot(X, FY)

limit 的值越大，图像的表示效果越好，但代价是运行ning 时间大得多。

在 matplotlib 中绘制积分奇异函数时出现的问题

Problems when plotting integrated singular functions in matplotlib

python

math

numpy

matplotlib

scientific-computing