Python: numpy.linalg.linalg.LinAlgError: Matrice singulière

Question

J'essaie de résoudre un système d'équations différentielles en utilisant ode solver dans scipy.integrate. Mon problème est, je reçois une erreur 'matrice singulière' quand ma matrice n'est pas supposée être singulière. Le problème principal est quand j'essaie de trouver l'inverse de la matrice B dans mon code ci-dessous. Dans le code ci-dessous B est une matrice 3x3, A est une matrice 3x1 et U est une matrice 3x1! Comment puis-je remédier à ce problème?

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

Pour les valeurs initiales avec beta=phi=0 et, en outre a[0]=0 de réglage a[k]=k ļ votre matrice B comporte une première ligne de zéros dans la première étape. Ce qui fait une matrice singulière.

Modification des données de test a[i]=i+1 supprime cette erreur, mais donne immédiatement une nouvelle erreur

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) 

qui est probablement parce que toutes vos dot variables sont encore des tableaux, de sorte que le vecteur de dérivés de retour est une liste de np.array de taille 3 chacun.

La première spéculation ici est que vous croyez à tort que * est le produit matriciel. Cela s'applique uniquement si vous construisez les objets en tant que np.matrix. Avec np.array, utilisez le dot pour les produits matriciels-vectoriels et np.linalg.solve au lieu du produit avec l'inverse.

Et en effet, qui est le cas, utilisez les lignes suivantes au lieu

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

Après que vous trouvez que vous utilisez le nom de la variable r tant pour le tableau des composants et pour le solveur ODE. Réparer que vous obtenez toujours que la taille du pas devient trop petite.

---

(ajouter Jul/19) Utiliser selon les commentaires de la variante

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

j'ai changé la boucle d'intégration à

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

où le traitement automatisé des composants de la solution marques pour un meilleur code et un tracé est produit chaque fois que t passe une valeur entière. Le dernier graphique est

avec la singularité en alpha se produisant à t=3.292894674 avec des valeurs largement instables dans alphadot variant -1.2e06-+2.9e06 avec des variations de la taille de 1e-09 dans le temps.