Python：numpy.linalg.linalg.LinAlgError：特異行列

Question

私はscipy.integrateでodeソルバーを使って微分方程式の系を解こうとしています。私の問題は、私の行列が特異ではないと思われるときに「特異行列」というエラーが出てきたことです。主な問題は、私が下のコードで行列Bの逆行列を見つけようとしているときです。下のコードでは、Bは3x3行列、Aは3x1行列、Uは3x1行列です。 この問題を解決するにはどうすればよいですか？ beta=phi=0とさらにa[0]=0を使用して初期値について

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.integrate import odeint 
import math 
import parameter_projection as pp 
from scipy.integrate import ode 
import sympy as sm 
c=10 
k,k1,k2,k3=10,9,8,7 
eta=2.5 
gamma,gamma1,gamma2=2,3,10 
a=[] 
for i in range(30): 
    a.append(i) 
a=np.array(a) 

def aero(t,Y): 
    V,alpha,beta,p,q,r,phi,theta=Y[0],Y[1],Y[2],Y[3],Y[4],Y[5],Y[6],Y[7] 
    sg=np.cos(alpha)*np.cos(beta)*np.sin(theta)-np.sin(beta)*np.sin(phi)*np.cos(theta)-np.sin(alpha)*np.cos(beta)*np.cos(phi)*np.cos(theta) 
    smcg=np.cos(alpha)*np.sin(beta)*np.sin(theta)+np.cos(beta)*np.sin(phi)*np.cos(theta)-np.sin(alpha)*np.sin(beta)*np.cos(phi)*np.cos(theta) 
    cmcg=np.sin(theta)*np.sin(alpha)+np.cos(alpha)*np.cos(phi)*np.cos(theta) 
    #error 
    ev=V-np.sin(t) 
    ebeta=beta-np.sin(t) 
    etheta=theta-np.sin(t) 
    ethetad=q*np.cos(phi)-r*np.sin(phi)-np.cos(t) 
    sv,sbeta,stheta=ev,ebeta,etheta+ethetad 
    s=np.array([[sv],[sbeta],[stheta]]) 
    A=np.array([[-a[1]*V**2*np.sin(beta)-a[2]*V**2*np.sin(beta)-a[4]*np.sin(gamma)-np.cos(t)],[p*np.sin(alpha)-r*np.cos(alpha)+a[16]*V+a[15]*smcg-np.cos(t)],[ethetad+np.cos(phi)*a[10]*p*r+np.cos(phi)*a[6]*(r**2-p**2)+np.cos(phi)*a[20]*V**2-np.cos(phi)*a[21]*sg+-q*np.sin(phi)*p-q*np.sin(phi)*q*np.sin(phi)*np.tan(theta)-q*np.sin(phi)*r*np.cos(phi)*np.tan(theta)-np.sin(phi)*a[11]*p*q+a[12]*q*r-a[13]*V**2+r*np.cos(phi)*p+r*q*np.cos(phi)*np.sin(phi)*np.tan(theta)+(r*np.cos(phi))**2*np.tan(theta)-np.cos(t)]]) 
    B=np.array([[a[0]*np.cos(alpha)*np.sin(beta),a[7]*np.sin(beta),a[0]*np.sin(alpha)*np.cos(beta)],[-a[9]*np.cos(alpha)*np.sin(beta)/V,a[22]*np.cos(beta)/V,-a[9]*np.sin(alpha)*np.sin(beta)/V],[a[29]*np.cos(phi),a[26]*np.sin(alpha)*np.sin(beta)*np.cos(phi)-a[27]*np.sin(phi)*np.cos(alpha),-a[25]*np.cos(phi)]]) 
    C=np.linalg.inv(B)*A 
    U=(C-np.linalg.inv(B)*k*s-np.linalg.inv(B)*eta*np.tanh(s)) 
    Vdot=a[0]*np.cos(alpha)*np.sin(beta)*U[0]-a[1]*V**2*np.cos(beta)-a[2]*V**2*np.sin(beta)-a[3]*sg+a[7]*np.sin(beta)*U[1]+a[0]*np.sin(alpha)*np.cos(beta)*U[2] 
    alphadot=q-(p*np.cos(alpha)+r*np.sin(alpha))*np.sin(beta)/np.cos(beta)+a[4]*V-a[14]*cmcg-a[8]*np.sin(alpha)*U[0]/V+a[8]*np.cos(alpha)*U[2]/V 
    betadot=p*np.sin(alpha)-r*np.cos(alpha)+a[16]*V+a[17]*smcg-a[9]*np.cos(alpha)*np.sin(beta)*U[0]/V+a[22]*np.cos(beta)*U[1]/V-a[9]*np.sin(alpha)*np.sin(beta)*U[2]/V 
    pdot=a[5]*q*r+a[17]*p*q+a[18]*V**2-a[19]*smcg+a[23]*U[0]-a[28]*U[2]+a[24]*np.sin(alpha)*np.cos(beta)*U[1] 
    qdot=a[10]*p*r+a[6]*(r**2-p**2)+a[20]*V**2-a[21]*sg+a[29]*U[0]-a[25]*U[2]+a[26]*np.sin(alpha)*np.sin(beta)*U[1] 
    rdot=a[11]*p*q-a[12]*q*r+a[13]*V**2+a[27]*np.cos(alpha)*U[2] 
    phidot=p+q*np.sin(phi)*np.tan(theta)+r*np.cos(phi)*np.tan(theta) 
    thetadot=q*np.cos(phi)-r*np.sin(phi) 
    return[Vdot,alphadot,betadot,pdot,qdot,rdot,phidot,thetadot] 

y0=[0.01,0.2,0,0,0,0,0,0.1] 
t0=0 
V=[] 
alpha=[] 
beta=[] 
p=[] 
q=[] 
r=[] 
phi=[] 
theta=[] 
t=[] 
r=ode(aero).set_integrator('dopri5',method='bdf',nsteps=1e10) 
r.set_initial_value(y0,t0) 
t1=10 
dt=.01 
while r.successful() and r.t<t1: 
    r.integrate(r.t+dt) 
    V.append(r.y[0]) 
    alpha.append(r.y[1]) 
    beta.append(r.y[2]) 
    p.append(r.y[3]) 
    q.append(r.y[4]) 
    r.append(r.y[5]) 
    phi.append(r.y[6]) 
    theta.append(r.y[7]) 
    t.append(r.t) 
V=np.array(V) 
alpha=np.array(alpha) 
beta=np.array(beta) 
p=np.array(p) 
q=np.array(q) 
r=np.array(r) 
phi=np.array(phi) 
theta=np.array(theta) 
plt.plot(t,V) 
plt.show() 

Answer 1

あなたマトリックスB°a[k]=k設定から最初のステップにおいてゼロの最初のラインを有しています。これは特異行列を作る。

a[i]=i+1にテストデータを変更すると、このエラーが削除されますが、すぐにすべてのあなたのdot変数がまだ配列されているので、返されたデリバティブのベクトルは大きさのnp.arrayのリストであるため、おそらく、新しいエラー

rv_cb_arr is NULL 
Call-back cb_fcn_in___user__routines failed. 
Traceback (most recent call last): 
    File "odeint-LA-singular_so45165407.py", line 60, in <module> 
     r.integrate(r.t+dt) 

を与えますそれぞれ3。

ここでの最初の推測は、誤って*がマトリックス製品であると誤解していることです。これは、オブジェクトをnp.matrixとして構成する場合にのみ適用されます。 np.arrayを使用すると、マトリックスベクトルプロダクトにはdot、逆行列プロダクトにはnp.linalg.solveが使用されます。

そして実際そうであることは、あなたがコンポーネントアレイのとODEソルバのために両方の変数名rを使用することを見つけるその後、次の行の代わりに

C = np.linalg.solve(B,A) 
    U = C-np.linalg.solve(B,k*s - eta*np.tanh(s)) 

を使用しています。あなたはまだステップサイズが小さくなりすぎるということを修復します。変異体を、コメント通り

U = -C-np.linalg.solve(B,k*s - eta*np.tanh(s)) 

を使用

---

（7月/ 19を追加） Iは

ソリューションコンポーネントの自動処理を行う

t0=0. 
y0=[0.01,0.2,0.,0.,0.,0.,0.,0.1] 
names = [ "V", "alpha", "beta", "p", "q", "r", "phi", "theta" ] 

sol = [ [yy] for yy in y0 ] 
t=[t0] 
solver=ode(aero).set_integrator('dopri5',nsteps=1e3) 
solver.set_initial_value(y0,t0) 
t1=3 
dt=5e-3 
for t1 in range(1,11): 
    if not solver.successful(): break; 
    while solver.successful() and solver.t<t1: 
     solver.integrate(solver.t+dt) 
     for k in range(8): sol[k].append(solver.y[k]); 
     t.append(solver.t) 
     print "step" 
    print len(t) 

    for k in range(8): 
     plt.subplot(4,2,k+1) 
     plt.plot(t,sol[k]) 
     plt.ylabel(names[k]) 
    plt.show() 

に統合ループを変更tが整数値を渡すたびにプロットが生成されます。最終的なグラフはalphaにおける特異点は、時間で1e-09の大きさの変化に-1.2e06から+2.9e06に変化alphadotで広く不安定な値とt=3.292894674で起こると

あります。