質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

2回答

2103閲覧

【Python】行列式=0となるような代入値の求め方

Masa06

総合スコア2

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2020/08/01 12:43

前提・実現したいこと

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])

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

can110

2020/08/01 13:54

omega=0じゃダメですか?
guest

回答2

0

ベストアンサー

以前に私が回答したこの質問と全く同じなので、参考にしてください。
【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

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

0

行列内に変数を含む

代数計算をしたいということでしたら、SymPyを使うのが良かろうかと思います。
Python (SymPy) による行列の計算

投稿2020/08/01 13:43

jeanbiego

総合スコア3966

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問