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

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

ただいまの
回答率

89.07%

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

解決済

回答 2

投稿

  • 評価
  • クリップ 0
  • VIEW 129

Masa06

score 1

前提・実現したいこと

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

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

作成中のコード

import cmath
import numpy as np
import scipy 
from scipy.integrate import quad 

#数値は仮
r0 = 1
r1 = 1
r2 = 1
r3 = 1
vr2 = 1
beta2 = 1
n = 1
Q = 1

ganma = - Q / cmath.tan(beta2)

#積分1
def M(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

A_M = M(lambda r: r2 / (r * cmath.sin(beta2)), r2, r1)

#複素積分計算1
def F_r2(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

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) 

#複素積分計算2
def F_r3(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

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) 

#複素積分計算3
def G_r2(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

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) 

#複素積分計算4
def G_r3(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

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) 

#raund_theta算出
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]

#各種値計算
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.matrix([
[a11, a12, a13, a14],
[0, a22, a23, a24],
[a31, a32, 0, 0],
[0, 0, a43, a44]
])
  • 気になる質問をクリップする

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

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

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • can110

    2020/08/01 22:54

    omega=0じゃダメですか?

    キャンセル

回答 2

checkベストアンサー

+1

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

複素数でとく場合

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

import cmath
import numpy as np
import scipy 
from scipy.integrate import quad 
import numpy.linalg as LA
from scipy import optimize

#積分1
def M(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])


#複素積分計算1
def F_r2(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])


#複素積分計算2
def F_r3(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])


#複素積分計算3
def G_r2(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])


#複素積分計算4
def G_r3(func, a, b, **kwargs):
    def real_func(r):
        return scipy.real(func(r))
    def imag_func(r):
        return scipy.imag(func(r))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

def A(x):
    #数値は仮
    r0 = 1
    r1 = 1
    r2 = 1
    r3 = 1
    vr2 = 1
    beta2 = 1
    n = 1
    Q = 1
    ganma = - Q / cmath.tan(beta2)

    omega = x[0]+x[1]*1j

    A_M = M(lambda r: r2 / (r * cmath.sin(beta2)), r2, r1)
    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) 
    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) 
    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) 
    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) 

    #raund_theta算出
    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]
    #各種値計算
    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.matrix([
    [a11, a12, a13, a14],
    [  0, a22, a23, a24],
    [a31, a32,   0,   0],
    [  0,   0, a43, a44]
    ])
    return np.abs(LA.det(A))

optimize.minimize(A, [0,1], method='Nelder-Mead')

結果

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

 final_simplex: (array([[-1.00106837e-05, -1.67414895e-05],
       [-1.63706702e-05,  4.94895503e-05],
       [-5.72935526e-05, -4.18866985e-06]]), array([1.47572037e-05, 3.94333956e-05, 4.34533398e-05]))
           fun: 1.475720371432839e-05
       message: 'Optimization terminated successfully.'
          nfev: 76
           nit: 41
        status: 0
       success: True
             x: array([-1.00106837e-05, -1.67414895e-05])

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/08/02 01:21

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

    キャンセル

  • 2020/08/02 02:47

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

    キャンセル

  • 2020/08/02 11:30

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

    キャンセル

0

行列内に変数を含む

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

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

  • ただいまの回答率 89.07%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る