行列の計算結果をグラフにしたい

前提・実現したいこと

S11とS21の計算結果を、グラフしたいですしたいです。S11とS21の計算過程で行列を用いています。今、fが1~10[GHz]の間でプロットしたいです。fが1つの値に決めるときは、計算ができるのですが、fをnp.linspace()などで表すとエラーが出ます。電磁波の計算をしています。F行列を求めてS行列を計算しています。

発生している問題・エラーメッセージ

TypeError: only length-1 arrays can be converted to Python scalars

エラーメッセージ

該当のソースコード

import numpy as np
import math
import matplotlib.pyplot as plt

def Q(f):
    return np.matrix([[1,Z/2],
                      [0,1]])
def W(f):
    return np.matrix([[np.cos((k*d)/2),Z0*np.sin((k*d)/2)*1j],
                      [Y0*np.sin((k*d)/2)*1j,np.cos((k*d)/2)]])

def E(f):
    return np.matrix([[1,0],
                      [Y,1]])

def FF(f):
    return Q(f)@W(f)@E(f)@W(f)@Q(f)

def F(f):
    n=5
    return np.linalg.matrix_power(FF, n)

f_min=1000000000
f_max=10000000000
f_points=100

f=np.linspace(f_min,f_max,f_points)

w=2*np.pi*f
d=5*10**(-3)
Z0=78
Y0=1/Z0
f0=2.4*10**9
Wse=2*np.pi*f0
Wsh=2*np.pi*f0
Cl=1.5*10**(-12)
Lr=1/(Wse**2*Cl)
Cr=Lr/(Z0**2)
Ll_=1/(Wsh**2*Cr)
e0=8.8541878128*(10**(-12))
u0=1.25663706212*(10**(-6))
k=2*np.pi*f*np.sqrt((Cr/d)*(Lr/d))
theta=(k*d)/2
w0=2*np.pi*f0
eeff=2.047
b0=w0*np.sqrt(u0*e0*eeff)
#ωはω0(共振周波数)のとき
l=1/b0*np.arctan(w0*Ll_/Z0)

b=w*np.sqrt(u0*e0*eeff)
Ll=Z0*np.tan(b*l)/w

Z=-1/(w*Cl)*1j
Y=-1/(w*Ll)*1j

matQ=np.zeros((f_points,2,2),dtype=complex)
matW=np.zeros((f_points,2,2),dtype=complex)
matE=np.zeros((f_points,2,2),dtype=complex)
matFF=np.zeros((f_points,2,2),dtype=complex)
matF=np.zeros((f_points,2,2),dtype=complex)

s11=np.zeros((f_points),dtype=complex)
s22=np.zeros((f_points),dtype=complex)

for i in range(f_points):
    matQ[i]=Q(f[i])
    matW[i]=W(f[i])
    matE[i]=E(f[i])
    matFF[i]=FF(f[i])
    matF[i]=F(f[i])
    
    A[i]=F[i[0,0]]
    B[i]=F[i[0,1]]
    C[i]=F[i[1,0]]
    D[i]=F[i[1,1]]
    
    s11[i]=(A[i]+B[i]/Z0-(Z0*C[i]+D[i]))/(A[i]+B[i]/Z0+Z0*C[i]+D[i])
    s21[i]=2/(A[i]+B[i]/Z0+Z0*C[i]+D[i])
    
S11=20*np.log10(abs(s11))
S21=20*np.log10(abs(s21))

plt.plot(f,S11,label="S11")
plt.plot(f,S21,label="S21")


plt.grid(True)
plt.xlabel("f",fontsize=20)
plt.ylabel("S-para",fontsize=20)
plt.title("S-para")
plt.tight_layout()
plt.show()ソースコード

試したこと

import numpy as np

f=2.010**9
w=2np.pi*f

d=510**(-3)
Z0=78
Y0=1/Z0
f0=2.4109
Wse=2np.pif0
Wsh=2np.pif0
Cl=1.5*10(-12)
Lr=1/(Wse2*Cl)
Cr=Lr/(Z02)
Ll_=1/(Wsh2Cr)
e0=8.8541878128(10(-12))
u0=1.25663706212*(10**(-6))
k=2np.pifnp.sqrt((Cr/d)(Lr/d))
theta=(kd)/2
w0=2np.pif0
eeff=2.047
b0=w0np.sqrt(u0e0eeff)
l=1/b0np.arctan(w0Ll_/Z0)

b=wnp.sqrt(u0e0eeff)
Ll=Z0np.tan(b*l)/w

Z=-1/(w*Cl)1j
Y=-1/(wLl)*1j

Q=np.matrix([[1,Z/2],
[0,1]])
W=np.matrix([[np.cos((kd)/2),Z0np.sin((kd)/2)1j],
[Y0np.sin((kd)/2)1j,np.cos((kd)/2)]])
E=np.matrix([[1,0],
[Y,1]])

FF=Q@W@E@W@Q

FF_=FF[0,0]*FF[1,1]-FF[0,1]*FF[1,0]
print("FF_={:.5f}".format(FF_))

n=5
F=np.linalg.matrix_power(FF, n)

A=F[0,0]
B=F[0,1]
C=F[1,0]
D=F[1,1]

s11=(A+B/Z0-(Z0C+D))/(A+B/Z0+Z0C+D)
S11=20*np.log10(abs(s11))

s21=2/(A+B/Z0+Z0C+D)
S21=20np.log10(abs(s21))

print("s11={:.5f}".format(s11))
print("S11={:.5f}".format(S11))
print("s21={:.5f}".format(s21))
print("S21={:.5f}".format(S21))

上記のように、fが1つの値を持つときは、計算できるのですが、S11とS21のグラフを作ろうとするとうまくいきません。横軸f、縦軸S11,S21で考えたいです。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

melian

2021/12/16 07:16

こちらの環境(Python 3.9.5/Numpy 1.21.4)で実行してみると以下のエラーが発生します。 matQ[i]= Q(f[i]) TypeError: only length-1 arrays can be converted to Python scalars

jbpb0

2021/12/16 08:08

> A[i]=F[i[0,0]] B[i]=F[i[0,1]] C[i]=F[i[1,0]] D[i]=F[i[1,1]] の「i[0,0]」とかって、何を意図してるのでしょうか？「i」って、「for i in range(f_points):」で作られる 0〜f_points-1 の整数の内のどれか一つだけが入ってる変数で、配列ではないので、 print(i[0,0]) とやったら、「TypeError: 'int' object is not subscriptable」というエラー出ますよ

kyoukyou

2021/12/16 08:10

ありがとうございます。申し訳ありません。自分のミスです。僕も、今確認すると同じエラーが発生しました。アドバイスなどがあればよろしくお願いします。

kyoukyou

2021/12/16 08:20

ありがとうございます。 A[i]=F[i[0,0]] B[i]=F[i[0,1]] C[i]=F[i[1,0]] D[i]=F[i[1,1]] A,B,C,Dは行列の要素です。i番目の行列Fの要素を取り出したいと考えてこのように書きました。ご指摘ありがとうございます。参考になります。

jbpb0

2021/12/16 08:22

あと、「def Q(f):」は、「f」とは無関係にいつも同じ結果が返るのですが、それは意図通りなのでしょうか？関数内で計算に使ってる「Z」には「f」と同じ数の要素があるので、もしかしたら「for i in range(f_points):」で作られる「i」番目の要素だけ使って計算した結果を返したい、のではないのでしょうか？「def W(f):」や「def E(f):」も同様です

jbpb0

2021/12/16 08:28

> i番目の行列Fの要素「F」は行列ではなくて関数です「def F(f):」ですから

kyoukyou

2021/12/16 08:35

その通りです。「i」番目同士の要素だけを用いて、S11,S21を計算して、その結果をプロットして曲線を描きたいと考えています。説明不足で申し訳ございません。

kyoukyou

2021/12/16 08:48

分かりました。 def F(f): は、使わない方がいいんですね。fをnp.linspaceで表して、それをx軸にして、S11,S21のグラフにしたいと考えました。なんどもありがとうございます。 import numpy as np import matplotlib.pyplot as plt f=np.linspace(1*10**9,10*10**10,100) w=2*np.pi*f d=5*10**(-3) Z0=78 Y0=1/Z0 f0=2.4*10**9 Wse=2*np.pi*f0 Wsh=2*np.pi*f0 Cl=1.5*10**(-12) Lr=1/(Wse**2*Cl) Cr=Lr/(Z0**2) Ll=1/(Wsh**2*Cr) e0=8.8541878128*(10**(-12)) u0=1.25663706212*(10**(-6)) k=2*np.pi*f*np.sqrt((Cr/d)*(Lr/d)) theta=(k*d)/2 w0=2*np.pi*f0 eeff=2.047 b0=w0*np.sqrt(u0*e0*eeff) #ωはω0(共振周波数)のとき l=1/b0*np.arctan(w0*Ll/Z0) Z=-1/(w*Ll)*1j Y=-1/(w*Cl)*1j #単位セルのabcd行列 a=(np.cos(theta))**2-(np.sin(theta))**2+Z/2*Y*(np.cos(theta))**2+(Y*Z0+Z*Y0)*np.sin(theta)*np.cos(theta)*1j b=Z*(np.cos(theta))**2-(Z+Z0**2*Y)*(np.sin(theta))**2+Z**2*Y/4*np.cos(theta)+(2*Z0+Z0*Z*Y+Z**2/(2*Z0))*np.sin(theta)*np.cos(theta)*1j c=Y*(np.cos(theta))**2+2*Y0*np.sin(theta)*np.cos(theta)*1j d=(np.cos(theta))**2-(np.sin(theta))**2+Z/2*Y*(np.cos(theta))**2+(Y*Z0+Z*Y0)*np.sin(theta)*np.cos(theta)*1j FF=np.matrix([[a,b],[c,d]]) # n セル数 n=5 F=FF*FF*FF*FF*FF A=F[0,0] B=F[0,1] C=F[1,0] D=F[1,1] #Sパラ #反射 s11=(A+B/Z0-(Z0*C+D))/(A+B/Z0+Z0*C+D) S11=20*np.log10(abs(s11)) #透過 s21=2/(A+B/Z0+Z0*C+D) S21=20*np.log10(abs(s21)) plt.plot(f,S11) plt.plot(f,S21) plt.grid(True) plt.xlabel("βd",fontsize=20) plt.ylabel("f",fontsize=20) plt.title("dispersion relation") plt.tight_layout() plt.show() 最初このようにしました。 ValueError: matrix must be 2-dimensional このようなエラーがでました。こちらをベースに考えなおしたほうがいいですかね。

jbpb0

2021/12/16 08:56

> 「i」番目同士の要素だけを用いてそうならば、このようにしたら、「Z」の「i」番目の要素だけが計算に使われます def Q(f): ↓ def Q(i): return np.matrix([[1,Z/2], ↓ return np.matrix([[1,Z[i]/2], 呼び出すのは、こうします matQ[i]=Q(f[i]) ↓ matQ[i]=Q(i) これで、「Z」の「i」番目の要素で計算された結果が、「matQ」の「i」番目の要素に入ります他の関数も同様で、関数の引数を「i」にして、「Z」の様に「f」と同じ要素数を持つ配列に「[i]」を付けます関数を呼び出す時の引数は「i」にしますあと、「関数F()」内の「return np.linalg.matrix_power(FF, n)」の「FF」は、「FF(i)」としないといけないような

jbpb0

2021/12/16 09:12

> def F(f): は、使わない方がいいんですね。使ってもいいですよ F()が2x2の行列を返す関数で、その結果の要素一つを取り出すのなら、こうします A[i]=F[i[0,0]] ↓ A[i]=F(i)[0,0] これで、「i」番目の要素で計算されたF(i)の結果の2x2の行列の内の左上の要素が、A[i]に格納されます

jbpb0

2021/12/16 09:16 編集

あと、配列A, B, C, Dを予め用意しておかないと、「for i in range(f_points):」のループ内で代入できませんもう一つ、下記が間違い s22=np.zeros((f_points),dtype=complex) ↓ s21=np.zeros((f_points),dtype=complex)

kyoukyou

2021/12/16 09:24

何度も本当にありがとうございます。配列A,B.C.Dの用意の仕方は A=np.zeros((f_points),dtype=complex) B=np.zeros((f_points),dtype=complex) C=np.zeros((f_points),dtype=complex) D=np.zeros((f_points),dtype=complex) このようだとだめですか？ TypeError: 'numpy.ndarray' object is not callable このようなエラーが表示されます。

melian

2021/12/16 09:29

質問内容を確定していただけると助かります。

jbpb0

2021/12/16 09:32 編集

A〜Dの定義は大丈夫です > TypeError: 'numpy.ndarray' object is not callable def FF(i): return Q(i)@W(i)@E(i)@W(i)@Q(i) としてますか？ matQ(i)とかしてませんか？

kyoukyou

2021/12/16 09:39

ご指摘いただいた、 A[i]=F[i[0,0]] ↓ A[i]=F(i)[0,0] これを A[i]=matF[i][0,0] このようにすると、できました。貴重なお時間ありがとうございます。

jbpb0

2021/12/16 09:40

これまでのコメントでいろいろ書いたことを(同様と書いたのを含めて)全部直したら、当方の手元のコードはエラーは出なくなりましたなので、見落としなくちゃんと直したら、少なくともエラーは出なくなります (計算が正しいかどうかは分かりませんが)

jbpb0

2021/12/16 09:42

> A[i]=F(i)[0,0] これを A[i]=matF[i][0,0] このようにすると、できました。あれ？関数F()を下記のように直したら、大丈夫なはずなのですが def F(i): n=5 return np.linalg.matrix_power(FF(i), n)

kyoukyou

2021/12/16 09:50

すみません。教えていただいた、 A[i]=F(i)[0,0] でもできてました。 A[i]=matF[i][0,0]とは何が違うのでしょうか？どちらでもできていました。僕の細かいミスがあったように思われます。本当にありがとうございました。今回の投稿は、きちんと作成します。はじめて利用したので、あまり使い方を把握していないのですが。

jbpb0

2021/12/16 10:07 編集

> A[i]=F(i)[0,0] でもできてました。 A[i]=matF[i][0,0]とは何が違うのでしょうか？結果は同じです関数F()の呼び出しが1回で済むので、「A[i]=matF[i][0,0]」の方がいいですね失礼しました【追記】回答のコードを、なるべく「mat*」を使うように書き直しました

行動規範の内容に同意します

回答1件

ベストアンサー

下記のように変えたらエラーが出なくなり、グラフが表示されましたので、ご参考までに
ただし、計算結果が正しいのかは分かりません

python
1import numpy as np
2import math
3import matplotlib.pyplot as plt
4
5
6#def Q(f):
7    #return np.matrix([[1,Z/2],
8def Q(i):
9    return np.matrix([[1,Z[i]/2],
10                      [0,1]])
11
12#def W(f):
13    #return np.matrix([[np.cos((k*d)/2),Z0*np.sin((k*d)/2)*1j],
14                      #[Y0*np.sin((k*d)/2)*1j,np.cos((k*d)/2)]])
15def W(i):
16    return np.matrix([[np.cos((k[i]*d)/2),Z0*np.sin((k[i]*d)/2)*1j],
17                      [Y0*np.sin((k[i]*d)/2)*1j,np.cos((k[i]*d)/2)]])
18
19#def E(f):
20def E(i):
21    return np.matrix([[1,0],
22                      #[Y,1]])
23                      [Y[i],1]])
24
25#def FF(f):
26    #return Q(f)@W(f)@E(f)@W(f)@Q(f)
27def FF(i):
28    return matQ[i]@matW[i]@matE[i]@matW[i]@matQ[i]
29
30#def F(f):
31def F(i):
32    n=5
33    #return np.linalg.matrix_power(FF, n)
34    return np.linalg.matrix_power(matFF[i], n)
35
36
37f_min=1000000000
38f_max=10000000000
39f_points=100
40
41f=np.linspace(f_min,f_max,f_points)
42
43w=2*np.pi*f
44d=5*10**(-3)
45Z0=78
46Y0=1/Z0
47f0=2.4*10**9
48Wse=2*np.pi*f0
49Wsh=2*np.pi*f0
50Cl=1.5*10**(-12)
51Lr=1/(Wse**2*Cl)
52Cr=Lr/(Z0**2)
53Ll_=1/(Wsh**2*Cr)
54e0=8.8541878128*(10**(-12))
55u0=1.25663706212*(10**(-6))
56k=2*np.pi*f*np.sqrt((Cr/d)*(Lr/d))
57theta=(k*d)/2
58w0=2*np.pi*f0
59eeff=2.047
60b0=w0*np.sqrt(u0*e0*eeff)
61#ωはω0(共振周波数)のとき
62l=1/b0*np.arctan(w0*Ll_/Z0)
63
64b=w*np.sqrt(u0*e0*eeff)
65Ll=Z0*np.tan(b*l)/w
66
67Z=-1/(w*Cl)*1j
68Y=-1/(w*Ll)*1j
69
70
71matQ=np.zeros((f_points,2,2),dtype=complex)
72matW=np.zeros((f_points,2,2),dtype=complex)
73matE=np.zeros((f_points,2,2),dtype=complex)
74matFF=np.zeros((f_points,2,2),dtype=complex)
75matF=np.zeros((f_points,2,2),dtype=complex)
76
77# 追加
78A=np.zeros((f_points),dtype=complex)
79B=np.zeros((f_points),dtype=complex)
80C=np.zeros((f_points),dtype=complex)
81D=np.zeros((f_points),dtype=complex)
82
83s11=np.zeros((f_points),dtype=complex)
84#s22=np.zeros((f_points),dtype=complex)
85s21=np.zeros((f_points),dtype=complex)
86
87for i in range(f_points):
88    #matQ[i]=Q(f[i])
89    #matW[i]=W(f[i])
90    #matE[i]=E(f[i])
91    #matFF[i]=FF(f[i])
92    #matF[i]=F(f[i])
93    matQ[i]=Q(i)
94    matW[i]=W(i)
95    matE[i]=E(i)
96    matFF[i]=FF(i)
97    matF[i]=F(i)
98    #A[i]=F[i[0,0]]
99    #B[i]=F[i[0,1]]
100    #C[i]=F[i[1,0]]
101    #D[i]=F[i[1,1]]
102    A[i]=matF[i][0,0]
103    B[i]=matF[i][0,1]
104    C[i]=matF[i][1,0]
105    D[i]=matF[i][1,1]
106    s11[i]=(A[i]+B[i]/Z0-(Z0*C[i]+D[i]))/(A[i]+B[i]/Z0+Z0*C[i]+D[i])
107    s21[i]=2/(A[i]+B[i]/Z0+Z0*C[i]+D[i])
108
109S11=20*np.log10(abs(s11))
110S21=20*np.log10(abs(s21))
111
112
113plt.plot(f,S11,label="S11")
114plt.plot(f,S21,label="S21")
115
116plt.grid(True)
117plt.xlabel("f",fontsize=20)
118plt.ylabel("S-para",fontsize=20)
119plt.title("S-para")
120plt.tight_layout()
121plt.show()