x<30でF(x)=a みたいにF(x)が変化するような関数をプロットしたい2

Question

条件によって変化するF(x)をプロットしたいがエラーが出てしまう。

from scipy.special import kv
import matplotlib.pyplot as plt
from scipy.integrate import quad
import numpy as np
import math
from math import gamma

xs = np.linspace(0, 100, 10000)
f = lambda z: kv(5/3,z)
F = [quad(f,x,np.inf)[0]*x for x in xs]

a = gamma(1/3) 
G = [(4*math.pi/np.sqrt(3)/a)*(x/2)**(1/3) for x in xs]
H = [((math.pi/2)**(1/2))*(x**(1/2))*(math.exp(-x))for x in xs]


def A(x):
    if x <= 5.0*1e-3:
        return G 
    elif 5.0*1-3 < x < 30:
        return F
    elif 30 <= x:
        return H

y = [A(i) for i in xs]

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.grid()
ax.plot(xs,y)
ax.set_xlim(1e-3, 1e1)
ax.set_ylim(1e-3, 1e0)
ax.set_yscale('log')  
ax.set_xscale('log')  
ax.set_title('log')
ax.set_xlabel('x')
ax.set_ylabel('')
plt.show()

エラー
TypeError                                 Traceback (most recent call last)
TypeError: float() argument must be a string or a number, not 'list'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-14-37ff64188e32> in <module>
     28 ax = fig.add_subplot(1,1,1)
     29 ax.grid()
---> 30 ax.plot(xs,y)
     31 ax.set_xlim(1e-3, 5*1e1)
     32 ax.set_ylim(1e-3, 1e0)

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\axes\_axes.py in plot(self, scalex, scaley, data, *args, **kwargs)
   1665         lines = [*self._get_lines(*args, data=data, **kwargs)]
   1666         for line in lines:
-> 1667             self.add_line(line)
   1668         self.autoscale_view(scalex=scalex, scaley=scaley)
   1669         return lines

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\axes\_base.py in add_line(self, line)
   1900             line.set_clip_path(self.patch)
   1901 
-> 1902         self._update_line_limits(line)
   1903         if not line.get_label():
   1904             line.set_label('_line%d' % len(self.lines))

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\axes\_base.py in _update_line_limits(self, line)
   1922         Figures out the data limit of the given line, updating self.dataLim.
   1923         """
-> 1924         path = line.get_path()
   1925         if path.vertices.size == 0:
   1926             return

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\lines.py in get_path(self)
   1025         """
   1026         if self._invalidy or self._invalidx:
-> 1027             self.recache()
   1028         return self._path
   1029 

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\lines.py in recache(self, always)
    673         if always or self._invalidy:
    674             yconv = self.convert_yunits(self._yorig)
--> 675             y = _to_unmasked_float_array(yconv).ravel()
    676         else:
    677             y = self._y

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\cbook\__init__.py in _to_unmasked_float_array(x)
   1388         return np.ma.asarray(x, float).filled(np.nan)
   1389     else:
-> 1390         return np.asarray(x, float)
   1391 
   1392 

C:\ProgramData\Anaconda3\lib\site-packages\numpy\core\_asarray.py in asarray(a, dtype, order)
     83 
     84     """
---> 85     return array(a, dtype, copy=False, order=order)
     86 
     87 

ValueError: setting an array element with a sequence.

from scipy.special import kv
import matplotlib.pyplot as plt
from scipy.integrate import quad
import numpy as np
import math
from math import gamma

xs = np.linspace(0, 50, 10000)
f = lambda z: kv(5/3,z)
F = [quad(f,x,np.inf)[0]*x for x in xs]

a = gamma(1/3) 
G = [(4*math.pi/np.sqrt(3)/a)*(x/2)**(1/3) for x in xs]
H = [((math.pi/2)**(1/2))*(x**(1/2))*(math.exp(-x))for x in xs]


def A(x):
    if x <= 5.0*1e-3:
        return (4*math.pi/np.sqrt(3)/a)*(x/2)**(1/3) 
    elif 5.0*1-3 < x < 30:
        return quad(f,x,np.inf)[0]*x
    elif 30 <= x:
        return ((math.pi/2)**(1/2))*(x**(1/2))*(math.exp(-x))

y = [A(i) for i in xs]

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.grid()
ax.plot(xs,y)
ax.set_xlim(1e-3, 5*1e1)
ax.set_ylim(1e-3, 1e0)
ax.set_yscale('log')  
ax.set_xscale('log')  
ax.set_title('log')
ax.set_xlabel('x')
ax.set_ylabel('')
plt.show()

のようにすると、
途中できれたものが出てくる

Answer 1

エラーに関して

最初のエラーに関してですが、内包表記が必要となる部分はxsからyを作る部分のみです。他の部分は普通に関数として定義するようにして下さい。

f = lambda z: kv(5 / 3, z)
a = gamma(1 / 3)


def F(x):
    return quad(f, x, np.inf)[0] * x


def G(x):
    return (4 * math.pi / np.sqrt(3) / a) * (x / 2) ** (1 / 3)


def H(x):
    return ((math.pi / 2) ** (1 / 2)) * (x ** (1 / 2)) * (math.exp(-x))


def A(x):
    if x <= 5.0 * 1e-3:
        return G(x)
    elif 5.0 * 1 - 3 < x < 30:
        return F(x)
    elif 30 <= x:
        return H(x)

y = [A(i) for i in xs]

途中で切れたものが出てくる件に関して

考えられる点は2点あります。

1. A(x)の値がNaN
2. xの範囲が不十分

1に関しては計算式があっているかどうかを確認して下さい。2に関しては、対数スケールに対して

xs = np.linspace(0, 50, 10000)

としているため、刻み幅が約0.005(より少し大きい)となっています。そうするとx <= 5.0 * 1e-3は1点(x=0の場合のみ)しか得られないことになります。これを避けるためにはnp.logspace()を使って対数スケールに対して等幅となるようにxsを作成して下さい。例えば10の-3乗から10の2乗まで10000点で等幅とするには

xs = np.logspace(-3, 2, 10000)

とします。

Comments

[hataki7]

勉強になります。
主原因として範囲設定において、1e-3<x<30のeの部分が抜けていたことによるものでした。

Answer 2

elif 5.0*1-3 < x < 30:
↓
elif 5.0*1e-3 < x < 30: