前提・実現したいこと

Pythonで行列式=0となるような変数の代入値を求めるようなプログラムを作成しています．

下記作成中コードにて，「det A=0となるようなomegaの値を求める」コードを書きたいのですが，行列内に変数を含む場合のコーディングが分からないため，アドバイス頂けないでしょうか．

作成中のコード

Python
1import cmath
2import numpy as np
3import scipy 
4from scipy.integrate import quad 
5
6#数値は仮
7r0 = 1
8r1 = 1
9r2 = 1
10r3 = 1
11vr2 = 1
12beta2 = 1
13n = 1
14Q = 1
15
16ganma = - Q / cmath.tan(beta2)
17
18#積分1
19def M(func, a, b, **kwargs):
20    def real_func(r):
21        return scipy.real(func(r))
22    def imag_func(r):
23        return scipy.imag(func(r))
24    real_integral = quad(real_func, a, b, **kwargs)
25    imag_integral = quad(imag_func, a, b, **kwargs)
26    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
27
28A_M = M(lambda r: r2 / (r * cmath.sin(beta2)), r2, r1)
29
30#複素積分計算1
31def F_r2(func, a, b, **kwargs):
32    def real_func(r):
33        return scipy.real(func(r))
34    def imag_func(r):
35        return scipy.imag(func(r))
36    real_integral = quad(real_func, a, b, **kwargs)
37    imag_integral = quad(imag_func, a, b, **kwargs)
38    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
39
40A_F_r2 = F_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / r2)**(n + 1), r3, r2) 
41
42#複素積分計算2
43def F_r3(func, a, b, **kwargs):
44    def real_func(r):
45        return scipy.real(func(r))
46    def imag_func(r):
47        return scipy.imag(func(r))
48    real_integral = quad(real_func, a, b, **kwargs)
49    imag_integral = quad(imag_func, a, b, **kwargs)
50    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
51
52A_F_r3 = F_r3(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / r3)**(n + 1), r3, r3) 
53
54#複素積分計算3
55def G_r2(func, a, b, **kwargs):
56    def real_func(r):
57        return scipy.real(func(r))
58    def imag_func(r):
59        return scipy.imag(func(r))
60    real_integral = quad(real_func, a, b, **kwargs)
61    imag_integral = quad(imag_func, a, b, **kwargs)
62    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
63
64A_G_r2 = G_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r2 / r)**(n - 1), r2, r2) 
65
66#複素積分計算4
67def G_r3(func, a, b, **kwargs):
68    def real_func(r):
69        return scipy.real(func(r))
70    def imag_func(r):
71        return scipy.imag(func(r))
72    real_integral = quad(real_func, a, b, **kwargs)
73    imag_integral = quad(imag_func, a, b, **kwargs)
74    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
75
76A_G_r3 = G_r3(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r3 / r)**(n - 1), r3, r2) 
77
78#raund_theta算出
79raund_theta = cmath.exp(-1j *(omega * cmath.pi * (r3**2 - r2**2) / Q - n * ganma * cmath.log(r3 / r2) / Q)) - ((n - 1) / r3 ) * 0.5 * A_G_r3[0]
80
81#各種値計算
82a11 = r2 * vr2
83
84a12 = (1 - cmath.exp(-2* beta2 * 1j)) * ( 1j * r2 * omega / cmath.tan(beta2) - 0.5 * A_M[0] * n * omega )
85
86a13 = r1**n * omega
87
88a14 = omega / r1**n
89
90a22 = 1 - cmath.exp(-2* 1j * beta2)
91
92a23 = r1**(n -1)
93
94a24 = 1 / r1**(n + 1)
95
96a31 = (1j * omega - ganma /(2 * cmath.pi * r3**2) * 1j * n + Q/(2 * cmath.pi * r3**2)) * (0.5 * A_F_r3[0] - 0.5 * A_G_r3[0] - 1j * (r2 / r3)**(n + 1) * cmath.exp(-2 * 1j * beta2) * (-1j *((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0])))) + Q/(2 * cmath.pi * r3**2) * (raund_theta + 1j * (n + 1) * r2**(n + 1) / r3**(n + 2) * cmath.exp(-2 * 1j * beta2) * (-1j * ((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0]))))
97
98a32 = (1j * omega - ganma /(2 * cmath.pi * r3**2) * 1j * n + Q/(2 * cmath.pi * r3**2)) * (-1j * (r3 / r2)**(n - 1) - 1j * (r2 / r3)**(n + 1) * cmath.exp(-2 * 1j * beta2)) + (Q/(2 * cmath.pi * r3)) * (-1j * (n - 1) * r3**(n - 2) / r2**(n - 1) + 1j * (n + 1) * r2**(n + 1) / r3**(n + 2) * cmath.exp(-2 * 1j * beta2))
99
100a43 = omega * r0**(n - 1)
101
102a44 = -omega / r0**(n + 1)
103
104#行列計算
105A = np.matrix([
106[a11, a12, a13, a14],
107[0, a22, a23, a24],
108[a31, a32, 0, 0],
109[0, 0, a43, a44]
110])

8524ba23

2020/08/01 13:54

omega=0じゃダメですか？

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

回答2件

ベストアンサー

以前に私が回答したこの質問と全く同じなので、参考にしてください。
【python】numpyとscipyで変数を含んだ行列式を計算したい

複素数でとく場合

ある複素数をZとすると、Z=0⇔|Z|=0なので、今回は|Z|を最小化するという問題になります。
minimizeを使います。

python
1import cmath
2import numpy as np
3import scipy 
4from scipy.integrate import quad 
5import numpy.linalg as LA
6from scipy import optimize
7
8#積分1
9def M(func, a, b, **kwargs):
10    def real_func(r):
11        return scipy.real(func(r))
12    def imag_func(r):
13        return scipy.imag(func(r))
14    real_integral = quad(real_func, a, b, **kwargs)
15    imag_integral = quad(imag_func, a, b, **kwargs)
16    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
17
18
19#複素積分計算1
20def F_r2(func, a, b, **kwargs):
21    def real_func(r):
22        return scipy.real(func(r))
23    def imag_func(r):
24        return scipy.imag(func(r))
25    real_integral = quad(real_func, a, b, **kwargs)
26    imag_integral = quad(imag_func, a, b, **kwargs)
27    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
28
29
30#複素積分計算2
31def F_r3(func, a, b, **kwargs):
32    def real_func(r):
33        return scipy.real(func(r))
34    def imag_func(r):
35        return scipy.imag(func(r))
36    real_integral = quad(real_func, a, b, **kwargs)
37    imag_integral = quad(imag_func, a, b, **kwargs)
38    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
39
40
41#複素積分計算3
42def G_r2(func, a, b, **kwargs):
43    def real_func(r):
44        return scipy.real(func(r))
45    def imag_func(r):
46        return scipy.imag(func(r))
47    real_integral = quad(real_func, a, b, **kwargs)
48    imag_integral = quad(imag_func, a, b, **kwargs)
49    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
50
51
52#複素積分計算4
53def G_r3(func, a, b, **kwargs):
54    def real_func(r):
55        return scipy.real(func(r))
56    def imag_func(r):
57        return scipy.imag(func(r))
58    real_integral = quad(real_func, a, b, **kwargs)
59    imag_integral = quad(imag_func, a, b, **kwargs)
60    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
61
62def A(x):
63    #数値は仮
64    r0 = 1
65    r1 = 1
66    r2 = 1
67    r3 = 1
68    vr2 = 1
69    beta2 = 1
70    n = 1
71    Q = 1
72    ganma = - Q / cmath.tan(beta2)
73    
74    omega = x[0]+x[1]*1j
75    
76    A_M = M(lambda r: r2 / (r * cmath.sin(beta2)), r2, r1)
77    A_F_r2 = F_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / r2)**(n + 1), r3, r2) 
78    A_F_r3 = F_r3(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / r3)**(n + 1), r3, r3) 
79    A_G_r2 = G_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r2 / r)**(n - 1), r2, r2) 
80    A_G_r3 = G_r3(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r3 / r)**(n - 1), r3, r2) 
81
82    #raund_theta算出
83    raund_theta = cmath.exp(-1j *(omega * cmath.pi * (r3**2 - r2**2) / Q - n * ganma * cmath.log(r3 / r2) / Q)) - ((n - 1) / r3 ) * 0.5 * A_G_r3[0]
84    #各種値計算
85    a11 = r2 * vr2
86    a12 = (1 - cmath.exp(-2* beta2 * 1j)) * ( 1j * r2 * omega / cmath.tan(beta2) - 0.5 * A_M[0] * n * omega )
87    a13 = r1**n * omega
88    a14 = omega / r1**n
89    a22 = 1 - cmath.exp(-2* 1j * beta2)
90    a23 = r1**(n -1)
91    a24 = 1 / r1**(n + 1)
92    a31 = (1j * omega - ganma /(2 * cmath.pi * r3**2) * 1j * n + Q/(2 * cmath.pi * r3**2)) * (0.5 * A_F_r3[0] - 0.5 * A_G_r3[0] - 1j * (r2 / r3)**(n + 1) * cmath.exp(-2 * 1j * beta2) * (-1j *((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0])))) + Q/(2 * cmath.pi * r3**2) * (raund_theta + 1j * (n + 1) * r2**(n + 1) / r3**(n + 2) * cmath.exp(-2 * 1j * beta2) * (-1j * ((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0]))))
93    a32 = (1j * omega - ganma /(2 * cmath.pi * r3**2) * 1j * n + Q/(2 * cmath.pi * r3**2)) * (-1j * (r3 / r2)**(n - 1) - 1j * (r2 / r3)**(n + 1) * cmath.exp(-2 * 1j * beta2)) + (Q/(2 * cmath.pi * r3)) * (-1j * (n - 1) * r3**(n - 2) / r2**(n - 1) + 1j * (n + 1) * r2**(n + 1) / r3**(n + 2) * cmath.exp(-2 * 1j * beta2))
94    a43 = omega * r0**(n - 1)
95    a44 = -omega / r0**(n + 1)
96
97    #行列計算
98    A = np.matrix([
99    [a11, a12, a13, a14],
100    [  0, a22, a23, a24],
101    [a31, a32,   0,   0],
102    [  0,   0, a43, a44]
103    ])
104    return np.abs(LA.det(A))
105
106optimize.minimize(A, [0,1], method='Nelder-Mead')

結果

どういう計算をしているかわからないので、計算結果があっているかどうかはわかりません。
結果としてはomega = -1e-5 - 1.67e-5iとなります。

text
1 final_simplex: (array([[-1.00106837e-05, -1.67414895e-05],
2       [-1.63706702e-05,  4.94895503e-05],
3       [-5.72935526e-05, -4.18866985e-06]]), array([1.47572037e-05, 3.94333956e-05, 4.34533398e-05]))
4           fun: 1.475720371432839e-05
5       message: 'Optimization terminated successfully.'
6          nfev: 76
7           nit: 41
8        status: 0
9       success: True
10             x: array([-1.00106837e-05, -1.67414895e-05])

投稿2020/08/01 14:20

編集2020/08/01 16:44

Penpen7

総合スコア698

Masa06

2020/08/01 16:00

ご回答頂きありがとうございます．正に求めていたご回答でした．追加で質問させていただきたいのですがよろしいでしょうか．参考のURLを元にコードを書いたところ，「TypeError: must be real number, not function」とエラーが出てしまいました．エラーの内容から，実数解とするように言われていると察しましたが，複素数解として求めることは可能でしょうか？質問文のコード最下部を以下の通り修正したため，記載させていただきます．　※omegaの含まれるコードをdef()とすること，optimizeなどをimportするコードは記載済みです． ''' def A(omega): A = np.array([[a11, a12, a13, a14], [0, a22, a23, a24], [a31, a32, 0, 0], [0, 0, a43, a44]], dtype='complex64') return LA.det(A) print(optimize.fsolve(A, 0+0*1j)) '''

Penpen7

2020/08/01 16:14

複素数などは抜きにして、ひとまずAの中にもa11などを入れないといけませんよ。 def A(x): 以下インデントを入れる omega = x[0] a11 = r2 * vr2 a12 = (1 - cmath.exp(-2* beta2 * 1j)) * ( 1j * r2 * omega / cmath.tan(beta2) - 0.5 * A_M[0] * n * omega ) a13 = r1**n * omega a14 = omega / r1**n a22 = 1 - cmath.exp(-2* 1j * beta2) a23 = r1**(n -1) a24 = 1 / r1**(n + 1) a31 = (1j * omega - ganma /(2 * cmath.pi * r3**2) * 1j * n + Q/(2 * cmath.pi * r3**2)) * (0.5 * A_F_r3[0] - 0.5 * A_G_r3[0] - 1j * (r2 / r3)**(n + 1) * cmath.exp(-2 * 1j * beta2) * (-1j *((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0])))) + Q/(2 * cmath.pi * r3**2) * (raund_theta + 1j * (n + 1) * r2**(n + 1) / r3**(n + 2) * cmath.exp(-2 * 1j * beta2) * (-1j * ((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0])))) a32 = (1j * omega - ganma /(2 * cmath.pi * r3**2) * 1j * n + Q/(2 * cmath.pi * r3**2)) * (-1j * (r3 / r2)**(n - 1) - 1j * (r2 / r3)**(n + 1) * cmath.exp(-2 * 1j * beta2)) + (Q/(2 * cmath.pi * r3)) * (-1j * (n - 1) * r3**(n - 2) / r2**(n - 1) + 1j * (n + 1) * r2**(n + 1) / r3**(n + 2) * cmath.exp(-2 * 1j * beta2)) a43 = omega * r0**(n - 1) a44 = -omega / r0**(n + 1) A = np.array([[a11, a12, a13, a14], [0, a22, a23, a24], [a31, a32, 0, 0], [0, 0, a43, a44]], dtype='complex64') return LA.det(A)

Penpen7

2020/08/01 16:21

fsolveはおそらく実数でしか解けないと思います。ということは、実数omegaの実部と、虚部を別の文字として解釈して、np.abs(LA.det(A))が0になるようにとく必要があると思います。

Penpen7

2020/08/01 17:47

omega=0を代入してしまうと|A|=0となってしまうため、omegaがそちらの方向に向かっている気がします。

Masa06

2020/08/02 02:30

コードまで修正していただき，ご丁寧にご回答いただきありがとうございました．

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