私はscipy.integrateでodeソルバーを使って微分方程式の系を解こうとしています。私の問題は、私の行列が特異ではないと思われるときに「特異行列」というエラーが出てきたことです。主な問題は、私が下のコードで行列Bの逆行列を見つけようとしているときです。下のコードでは、Bは3x3行列、Aは3x1行列、Uは3x1行列です。 この問題を解決するにはどうすればよいですか? beta=phi=0
とさらにa[0]=0
を使用して初期値についてPython:numpy.linalg.linalg.LinAlgError:特異行列
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()
こんにちは、私はあなたのソリューションを試しましたが、私はrからsにソリューション名を変更しましたが、ステップサイズは非常に小さくなっています。私はまだグラフを得ることができません。 –
'dt = 1e-8'を設定し、その導関数が無限になることを見るために' alpha'をプロットします。つまり、 'alpha'は特異点が' t_sing = 4.50321260794985854e-06'。その行動を推進しているのはすぐには見えません。 – LutzL
こんにちは! LutzL実際にコードにエラーがありました。私はUの表現を-C-np.linalg.solve(B、ks-etanp.tanh(s))に変更しました。また、k *とη* np.tanh(s)の間にマイナスの代わりにプラス記号があると思います。 しかし、実際には過去5〜10分間実行しているので、コードを実行するのに時間がかかります! –