teratail
回答率
85.35%
teratail header banner
teratail header banner
質問するログイン新規登録
トップPythonに関する質問【Python】行列式=0となるような代入値の求め方

編集履歴

回答編集履歴

1

コード例追加

2020/08/01 16:44

投稿

Penpen7
Penpen7

スコア698

test CHANGED
@@ -1,3 +1,261 @@
1
1
  以前に私が回答したこの質問と全く同じなので、参考にしてください。
2
2
 
3
3
  [【python】numpyとscipyで変数を含んだ行列式を計算したい](https://teratail.com/questions/275887#reply-393270)
4
+
5
+
6
+
7
+ # 複素数でとく場合
8
+
9
+ ある複素数をZとすると、Z=0⇔|Z|=0なので、今回は|Z|を最小化するという問題になります。
10
+
11
+ ``minimize``を使います。
12
+
13
+
14
+
15
+ ```python
16
+
17
+ import cmath
18
+
19
+ import numpy as np
20
+
21
+ import scipy
22
+
23
+ from scipy.integrate import quad
24
+
25
+ import numpy.linalg as LA
26
+
27
+ from scipy import optimize
28
+
29
+
30
+
31
+ #積分1
32
+
33
+ def M(func, a, b, **kwargs):
34
+
35
+ def real_func(r):
36
+
37
+ return scipy.real(func(r))
38
+
39
+ def imag_func(r):
40
+
41
+ return scipy.imag(func(r))
42
+
43
+ real_integral = quad(real_func, a, b, **kwargs)
44
+
45
+ imag_integral = quad(imag_func, a, b, **kwargs)
46
+
47
+ return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
48
+
49
+
50
+
51
+
52
+
53
+ #複素積分計算1
54
+
55
+ def F_r2(func, a, b, **kwargs):
56
+
57
+ def real_func(r):
58
+
59
+ return scipy.real(func(r))
60
+
61
+ def imag_func(r):
62
+
63
+ return scipy.imag(func(r))
64
+
65
+ real_integral = quad(real_func, a, b, **kwargs)
66
+
67
+ imag_integral = quad(imag_func, a, b, **kwargs)
68
+
69
+ return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
70
+
71
+
72
+
73
+
74
+
75
+ #複素積分計算2
76
+
77
+ def F_r3(func, a, b, **kwargs):
78
+
79
+ def real_func(r):
80
+
81
+ return scipy.real(func(r))
82
+
83
+ def imag_func(r):
84
+
85
+ return scipy.imag(func(r))
86
+
87
+ real_integral = quad(real_func, a, b, **kwargs)
88
+
89
+ imag_integral = quad(imag_func, a, b, **kwargs)
90
+
91
+ return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
92
+
93
+
94
+
95
+
96
+
97
+ #複素積分計算3
98
+
99
+ def G_r2(func, a, b, **kwargs):
100
+
101
+ def real_func(r):
102
+
103
+ return scipy.real(func(r))
104
+
105
+ def imag_func(r):
106
+
107
+ return scipy.imag(func(r))
108
+
109
+ real_integral = quad(real_func, a, b, **kwargs)
110
+
111
+ imag_integral = quad(imag_func, a, b, **kwargs)
112
+
113
+ return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
114
+
115
+
116
+
117
+
118
+
119
+ #複素積分計算4
120
+
121
+ def G_r3(func, a, b, **kwargs):
122
+
123
+ def real_func(r):
124
+
125
+ return scipy.real(func(r))
126
+
127
+ def imag_func(r):
128
+
129
+ return scipy.imag(func(r))
130
+
131
+ real_integral = quad(real_func, a, b, **kwargs)
132
+
133
+ imag_integral = quad(imag_func, a, b, **kwargs)
134
+
135
+ return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
136
+
137
+
138
+
139
+ def A(x):
140
+
141
+ #数値は仮
142
+
143
+ r0 = 1
144
+
145
+ r1 = 1
146
+
147
+ r2 = 1
148
+
149
+ r3 = 1
150
+
151
+ vr2 = 1
152
+
153
+ beta2 = 1
154
+
155
+ n = 1
156
+
157
+ Q = 1
158
+
159
+ ganma = - Q / cmath.tan(beta2)
160
+
161
+
162
+
163
+ omega = x[0]+x[1]*1j
164
+
165
+
166
+
167
+ A_M = M(lambda r: r2 / (r * cmath.sin(beta2)), r2, r1)
168
+
169
+ 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)
170
+
171
+ 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)
172
+
173
+ 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)
174
+
175
+ 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)
176
+
177
+
178
+
179
+ #raund_theta算出
180
+
181
+ 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]
182
+
183
+ #各種値計算
184
+
185
+ a11 = r2 * vr2
186
+
187
+ a12 = (1 - cmath.exp(-2* beta2 * 1j)) * ( 1j * r2 * omega / cmath.tan(beta2) - 0.5 * A_M[0] * n * omega )
188
+
189
+ a13 = r1**n * omega
190
+
191
+ a14 = omega / r1**n
192
+
193
+ a22 = 1 - cmath.exp(-2* 1j * beta2)
194
+
195
+ a23 = r1**(n -1)
196
+
197
+ a24 = 1 / r1**(n + 1)
198
+
199
+ 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]))))
200
+
201
+ 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))
202
+
203
+ a43 = omega * r0**(n - 1)
204
+
205
+ a44 = -omega / r0**(n + 1)
206
+
207
+
208
+
209
+ #行列計算
210
+
211
+ A = np.matrix([
212
+
213
+ [a11, a12, a13, a14],
214
+
215
+ [ 0, a22, a23, a24],
216
+
217
+ [a31, a32, 0, 0],
218
+
219
+ [ 0, 0, a43, a44]
220
+
221
+ ])
222
+
223
+ return np.abs(LA.det(A))
224
+
225
+
226
+
227
+ optimize.minimize(A, [0,1], method='Nelder-Mead')
228
+
229
+ ```
230
+
231
+
232
+
233
+ # 結果
234
+
235
+ どういう計算をしているかわからないので、計算結果があっているかどうかはわかりません。
236
+
237
+ 結果としてはomega = -1e-5 - 1.67e-5iとなります。
238
+
239
+ ```text
240
+
241
+ final_simplex: (array([[-1.00106837e-05, -1.67414895e-05],
242
+
243
+ [-1.63706702e-05, 4.94895503e-05],
244
+
245
+ [-5.72935526e-05, -4.18866985e-06]]), array([1.47572037e-05, 3.94333956e-05, 4.34533398e-05]))
246
+
247
+ fun: 1.475720371432839e-05
248
+
249
+ message: 'Optimization terminated successfully.'
250
+
251
+ nfev: 76
252
+
253
+ nit: 41
254
+
255
+ status: 0
256
+
257
+ success: True
258
+
259
+ x: array([-1.00106837e-05, -1.67414895e-05])
260
+
261
+ ```