Inverser une matrice en calcul formel

[Legion]

Bonjour,

Je suis actuellement en stage et je fais face à un petit problème avec Xcas.

Je souhaite calculer les R0 d'un modèles de propagations de deux maladies infectieuses à 7 états et 2 souches infectieuses.
Pour cela, j'utilise les outils d'Xcas permettant de faire du calcul formel.

Mon script tourne en boucle à partir du moment où je tente d'inverser une matrice 8x7 appelée V avec la fonction inv(V).
(V ne contient que des symboles)

Je viens de découvrir Xcas et malheureusement, je n'ai rien trouvé d'autre pour inverser une matrice.

Pouvez-vous m'aider s'il vous plaît?

Merci

[parisse]

Pouvez-vous donner les details? Ca permettra de comprendre precisement ce que vous entendez par inverser une matrice 8x7 (je suppose que c'est resoudre un systeme?) et voir la complexite du probleme (s'il s'agit d'inverser une matrice 7x7 avec plein de variables, rien d'etonnant a ce que ca bloque avec l'algorithme du pivot de Gauss par exemple)

[Legion]

Pouvez-vous donner les details? Ca permettra de comprendre precisement ce que vous entendez par inverser une matrice 8x7 (je suppose que c'est resoudre un systeme?) et voir la complexite du probleme (s'il s'agit d'inverser une matrice 7x7 avec plein de variables, rien d'etonnant a ce que ca bloque avec l'algorithme du pivot de Gauss par exemple)

Oui, effectivement, désolé d'avoir était aussi flou.

Tout d’abord, voici mon script :

Code : Tout sélectionner

purge(t_S, t_I1, t_I2, t_I12, t_IPrime1, t_IPrime2, t_I2_1, t_I1_2,
      g_1,g_2,g_3,g_4,g_5,g_6,g_7,g_8,g_9,g_10,g_11,g_12,g_13,g_14,g_15,g_16,g_17,g_18,g_19,g_20,g_21,g_22,g_23,g_24,g_25,g_26,g_27,g_28,g_29,g_30,g_31,g_32,g_33,g_34,g_35,g_36,g_37,g_38,g_39,g_40,g_41,g_42,g_43,g_44,g_45,g_46,g_47,g_48,
      S, I1, I2, I12, IPrime1, IPrime2, I2_1, I1_2,
      q1_S, q1_I1, q1_I2, q1_I12, q1_IPrime1, q1_IPrime2, q1_I2_1, q1_I1_2,
      q2_S, q2_I1, q2_I2, q2_I12, q2_IPrime1, q2_IPrime2, q2_I2_1, q2_I1_2,
      q12_S, q12_I1, q12_I2, q12_I12, q12_IPrime1, q12_IPrime2, q12_I2_1, q12_I1_2,
      a1_S, a1_I1, a1_I2, a1_I12, a1_IPrime1, a1_IPrime2, a1_I2_1, a1_I1_2,
      a2_S, a2_I1, a2_I2, a2_I12, a2_IPrime1, a2_IPrime2, a2_I2_1, a2_I1_2,
      a12_S, a12_I1, a12_I2, a12_I12, a12_IPrime1, a12_IPrime2, a12_I2_1, a12_I1_2,
      t_S_S, t_S_I1, t_S_I2, t_S_I12, t_S_IPrime1, t_S_IPrime2, t_S_I2_1, t_S_I1_2,
      t_I1_S, t_I1_I1, t_I1_I2, t_I1_I12, t_I1_IPrime1, t_I1_IPrime2, t_I1_I2_1, t_I1_I1_2,
      t_I2_S, t_I2_I1, t_I2_I2, t_I2_I12, t_I2_IPrime1, t_I2_IPrime2, t_I2_I2_1, t_I2_I1_2,
      t_I12_S, t_I12_I1, t_I12_I2, t_I12_I12, t_I12_IPrime1, t_I12_IPrime2, t_I12_I2_1, t_I12_I1_2,
      t_IPrime1_S, t_IPrime1_I1, t_IPrime1_I2, t_IPrime1_I12, t_IPrime1_IPrime1, t_IPrime1_IPrime2, t_IPrime1_I2_1, t_IPrime1_I1_2,
      t_IPrime2_S, t_IPrime2_I1, t_IPrime2_I2, t_IPrime2_I12, t_IPrime2_IPrime1, t_IPrime2_IPrime2, t_IPrime2_I2_1, t_IPrime2_I1_2,
      t_I2_1_S, t_I2_1_I1, t_I2_1_I2, t_I2_1_I12, t_I2_1_IPrime1, t_I2_1_IPrime2, t_I2_1_I2_1, t_I2_1_I1_2,
      t_I1_2_S, t_I1_2_I1, t_I1_2_I2, t_I1_2_I12, t_I1_2_IPrime1, t_I1_2_IPrime2, t_I1_2_I2_1, t_I1_2_I1_2,
      d1_S, d1_I1, d1_I2, d1_I12, d1_IPrime1, d1_IPrime2, d1_I2_1, d1_I1_2,
      d2_S, d2_I1, d2_I2, d2_I12, d2_IPrime1, d2_IPrime2, d2_I2_1, d2_I1_2,
      d12_S, d12_I1, d12_I2, d12_I12, d12_IPrime1, d12_IPrime2, d12_I2_1, d12_I1_2,
      l1, l2, l12,
      b1,b2,b12, theta, gamma, G, A, TWH, Dec, L, B
);






F1:= t_S*q1_S +  t_I1*q1_I1 + t_I2* q1_I2 + t_I12*q1_I12 + t_IPrime1*q1_IPrime1 + t_IPrime2*q1_IPrime2 + t_I2_1*q1_I2_1 + t_I1_2*q1_I1_2;

F2:= t_S*q2_S +  t_I1*q2_I1 + t_I2* q2_I2 + t_I12*q2_I12 + t_IPrime1*q2_IPrime1 + t_IPrime2*q2_IPrime2 + t_I2_1*q2_I2_1 + t_I1_2*q2_I1_2;

F12:= t_S*q12_S +  t_I1*q12_I1 + t_I2* q12_I2 + t_I12*q12_I12 + t_IPrime1*q12_IPrime1 + t_IPrime2*q12_IPrime2 + t_I2_1*q12_I2_1 + t_I1_2*q12_I1_2;


Psi_1 := (t_S * g_1 * t_I1 ) * S * B_1 * F1 * a1_S;
Psi_2 := (t_I1 * g_2 * t_I2 ) * I1 * d1_I1;
Psi_3 := (t_I1 * g_3 * t_I2 ) * I1 * B_2 * F2 * a2_I1;
Psi_4 := (t_I1 * g_4 * t_I2 ) * I1 * t_I1_I2;
Psi_5 := (t_I2 * g_5 * t_S ) * I2 * d2_I2;
Psi_6 := (t_S * g_6 * t_I2 ) * S * B_2 * F2 * a2_S;
Psi_7 := (t_I2 * g_7 * t_I1 ) * I2 * d2_I2;
Psi_8 := (t_I2 * g_8 * t_I1 ) * I2 * B_1 * F1 * a1_I2;
Psi_9 := (t_I2 * g_9 * t_I1 ) * I2 * t_I2_I1;
Psi_10 := (t_I1 * g_10 * t_S ) * I1 * d1_I1;
Psi_11 := (t_I1 * g_11 * t_I12 ) * I1 * B_2 * F2 * a2_I1;
Psi_12 := (t_S * g_12 * t_I12 ) * S * B_12 * F12 * a12_S;
Psi_13 := (t_I2 * g_13 * t_I12 ) * I2 * B_1 * F1 * a2_I2;   
Psi_14 := (t_I12 * g_14 * t_I1 ) * I12 * d2_I12;
Psi_15 := (t_I12 * g_15 * t_S ) * I12 * d12_I12;
Psi_16 := (t_I12 * g_16 * t_I2 ) * I12 * d1_I12;
Psi_17 := (t_I12 * g_17 * t_IPrime1 ) * I12 * d2_I12;
Psi_18 := (t_IPrime1 * g_18 * t_S ) * IPrime1 * d1_IPrime1;
Psi_19 := (t_I12 * g_19 * t_IPrime2 ) * I12 * d1_I12;
Psi_20 := (t_IPrime2 * g_20 * t_S ) * IPrime2 * d2_IPrime2;
Psi_21 := (t_S * g_21 * t_I1_2 ) * S * B_12 * F12 * a12_S;
Psi_22 := (t_I1 * g_22 * t_I1_2 ) * I1 * B_2 * F2 * a2_I1_2;
Psi_23 := (t_S * g_23 * t_I2_1 ) * S * B_12 * F12 * a12_S;
Psi_24 := (t_I2 * g_24 * t_I2_1 ) * I2 * B_2 * F2 * a1_I2;
Psi_25 := (t_I1_2 * g_25 * t_S ) * I1_2 * d12_I1_2;
Psi_26 := (t_I1_2 * g_26 * t_I1 ) * I1_2 * d2_I1_2;
Psi_27 := (t_I1_2 * g_27 * t_I2 ) * I1_2 * d1_I1_2;
Psi_28 := (t_I2_1 * g_28 * t_S ) * I2_1 * d12_I2_1;
Psi_29 := (t_I2_1 * g_29* t_I2 ) * I2_1 * d1_I2_1;
Psi_30 := (t_I2_1 * g_30 * t_I1 ) * I2_1 * d2_I2_1;;
Psi_31 := (t_I1_2 * g_31 * t_I12 ) * I1_2 * d1_I1_2;
Psi_32 := (t_I1_2 * g_32 * t_I12 ) * I12 * B_2 * F2 * a2_I1_2;
Psi_33 := (t_I1_2 * g_33 * t_I12 ) * I1 * t_I1_2_I12;
Psi_34 := (t_I12 * g_31 * t_I2_1 ) * I12 * d1_I12;
Psi_35 := (t_I12 * g_35 * t_I2_1 ) * I12 * B_2 * F2 * a2_I12;
Psi_36 := (t_I12 * g_36 * t_I2_1 ) * I12 * t_I12_I2_1;
Psi_37 := (t_I1_2 * g_37 * t_I2_2 ) * I1_2 * d1_I1_2;
Psi_38 := (t_I1_2 * g_38 * t_I2_1 ) * I1_2 * B_2 * F2 * a2_I1_2;
Psi_39 := (t_I1_2 * g_39 * t_I2_1 ) * I1_2 * t_I1_2_I2_1;
Psi_40 := (t_I2_1 * g_40 * t_I12 ) * I2_1 * d2_I2_1;
Psi_41 := (t_I2_1 * g_41 * t_I12 ) * I2_1 * B_1 * F1 * a1_I2_1;
Psi_42 := (t_I2_1 * g_42 * t_I12 ) * I2_1 * t_I2_1_I12;
Psi_43 := (t_I12 * g_40 * t_I1_2 ) * I12 * d2_I12;
Psi_44 := (t_I12 * g_44 * t_I1_2 ) * I12 * B_1 * F1 * a1_I12;
Psi_45 := (t_I12 * g_45 * t_I1_2 ) * I12 * t_I12_I1_2;
Psi_46 := (t_I2_1 * g_46 * t_I1_2 ) * I2_1 * d2_I2_1;
Psi_47 := (t_I2_1 * g_47 * t_I1_2 ) * I2_1 * B_1 * F1 * a1_I2_1;
Psi_48 := (t_I2_1 * g_48 * t_I1_2 ) * I2_1 * t_I2_1_I1_2;

afficher("Psi");

Vbase:=[- (Psi_7 +  Psi_9 + Psi_14 + Psi_26 + Psi_30 - Psi_2 - Psi_3 - Psi_4 - Psi_10 - Psi_11 - Psi_22 ),
                - (Psi_2 + Psi_4 +Psi_16 + Psi_27 + Psi_29 - Psi_5 - Psi_7 - Psi_8 - Psi_9 - Psi_13 - Psi_24), 
                - (Psi_31 + Psi_33 + Psi_40 + Psi_42 - Psi_14 - Psi_15 - Psi_16 - Psi_17 - Psi_19 - Psi_34 - Psi_35 - Psi_36 - Psi_43 - Psi_44 - Psi_45), 
                Psi_18, 
                Psi_20, 
                -(Psi_34 +  Psi_36 + Psi_37 + Psi_39 - Psi_28 - Psi_29 - Psi_30 - Psi_40 - Psi_41 - Psi_42 - Psi_46 - Psi_47 - Psi_48), 
                -(Psi_43 +Psi_45 + Psi_46 + Psi_48 - Psi_25 - Psi_26 - Psi_27 - Psi_31 - Psi_32 - Psi_33 - Psi_37 - Psi_38 - Psi_39)];

afficher("VBase");

V:=[[diff(Vbase[1],I1),diff(Vbase[1],I2),diff(Vbase[1],I12),diff(Vbase[1],IPrime1),diff(Vbase[1],IPrime2),diff(Vbase[1],I2_1),diff(Vbase[1],I1_2)],
         [diff(Vbase[2],I1),diff(Vbase[2],I2),diff(Vbase[2],I12),diff(Vbase[2],IPrime1),diff(Vbase[2],IPrime2),diff(Vbase[2],I2_1),diff(Vbase[2],I1_2)],
         [diff(Vbase[3],I1),diff(Vbase[3],I2),diff(Vbase[3],I12),diff(Vbase[3],IPrime1),diff(Vbase[3],IPrime2),diff(Vbase[3],I2_1),diff(Vbase[3],I1_2)],
         [diff(Vbase[4],I1),diff(Vbase[4],I2),diff(Vbase[4],I12),diff(Vbase[4],IPrime1),diff(Vbase[4],IPrime2),diff(Vbase[4],I2_1),diff(Vbase[4],I1_2)],
         [diff(Vbase[5],I1),diff(Vbase[5],I2),diff(Vbase[5],I12),diff(Vbase[5],IPrime1),diff(Vbase[5],IPrime2),diff(Vbase[5],I2_1),diff(Vbase[5],I1_2)],
         [diff(Vbase[6],I1),diff(Vbase[6],I2),diff(Vbase[6],I12),diff(Vbase[6],IPrime1),diff(Vbase[6],IPrime2),diff(Vbase[6],I2_1),diff(Vbase[6],I1_2)],
         [diff(Vbase[7],I1),diff(Vbase[7],I2),diff(Vbase[7],I12),diff(Vbase[7],IPrime1),diff(Vbase[7],IPrime2),diff(Vbase[7],I2_1),diff(Vbase[7],I1_2)]];

afficher("V");
invV:= inverse(V);

afficher("invV");
R:=[[r,0,0,0,0,0,0],
        [0,r,0,0,0,0,0],
        [0,0,r,0,0,0,0],
        [0,0,0,r,0,0,0],
        [0,0,0,0,r,0,0],
        [0,0,0,0,0,r,0],
        [0,0,0,0,0,0,r]];

afficher("R");
res := simplifier(developper(solve(det(F*invV-R)=0,r)));

afficher(res);

La fonction "afficher" sert uniquement à suivre l'avancée des calculs.
Le script s'arrête entre afficher("V") et afficher("invV").
J'en déduis donc que le problème vient de la ligne invV:= inverse(V).
J'avoue ne pas comprendre pourquoi après une journée de calculs, je n'ai pas eu de résultat.

Merci de votre aide.

[frederic han]

bonjour quelques conseils vu le grand nombre de variables.

il vaut mieux essayer petit a petit.

par exemple j'ai peut etre eu un pb de copier/coller, mais je trouve que le determiant de V est nul. vu qu'il y a quelques 0 et beaucoup de variables, la methode des mineurs est plus efficace:

Code : Tout sélectionner

det_minor(V)

ensuite, toujours pour eviter les denominateurs trop grand, il vaut peut etre mieux trouver la comatrice de V
j'ai reussi a faire:

Code : Tout sélectionner

L:=adjoint_matrix(V)

en 1min mais 730Mo
la comatrice de V est donc:

Code : Tout sélectionner

A:=L[1][6]

NB: si vous voulez un polynome caracteristique il vaut mieux le demander directement avec pcar.

[parisse]

Trop long pour inv, par contre on peut faire det(V-r*F) sans changer le probleme il me semble, ou plus efficacement det_minor(V-r*F) ce qui prend quelques minutes, mais ensuite ca m'etonnerait que solve puisse faire quoi que ce soit.
pour fred: interessant de voir que adjoint_matrix va plus vite!

[Legion]

Merci à tous pour votre aide.
Je vais essayer de faire det_minor(V-R*F) mais si solve ne peut rien faire alors je me demande si je ne suis pas tout simplement fichue ?

Pourriez vous me dire comment arrivez vous à trouver que "det(V-r*F) = det(F*invV-R)", s'il vous plaît ?

Comme je travaille sur un cluster, j'ai recherché s'il était possible de lancer un script Xcas en terminal sans avoir a passé par l'interface graphique. Malheureusement, je n'ai rien trouvé.
Avez vous une idée?

Merci

[parisse]

det(F*invV-R)=det((F-R*V)*invV)=det(F-r*V)*det(invV) donc les annulations des 2 sont equivalentes.
Pour executer en mode script, il suffit d'utiliser icas au lieu de xcas
icas nom_de_fichier
ou nom_de_fichier est le nom du fichier contenant vos commandes.

[Legion]

det(F*invV-R)=det((F-R*V)*invV)=det(F-r*V)*det(invV) donc les annulations des 2 sont equivalentes.
Pour executer en mode script, il suffit d'utiliser icas au lieu de xcas
icas nom_de_fichier
ou nom_de_fichier est le nom du fichier contenant vos commandes.

Merci.

J'ai modifié mon script en ajoutant :

Code : Tout sélectionner

deter := det_minor(V-R*F) ;
afficher("Deter");
res := solve(deter = 0,r);

Je confirme que solve est toujours en cours. Y aurait il un fonction plus efficace?

[parisse]

Non, et meme si la factorisation aboutit, ca donnerait surement des resultats beaucoup trop gros donc inutilisables au pire et tres inefficaces au mieux. Je pense que l'approche vouloir une solution analytique n'est pas la bonne (un peu comme vouloir une formule pour les equations du 4eme degre). Il serait sans doute beaucoup plus efficace de remplacer le script actuel par une fonction prenant en argument les valeurs des parametres et renvoyant le resultat du solve(...) dans ce cas.
Donc remplacer le mot clef purge par f puis mettre une accolade { apres la parenthese fermante a la place du point virgule et terminer par return res; }

[frederic han]

Merci à tous pour votre aide.

Pourriez vous me dire comment arrivez vous à trouver que "det(V-r*F) = det(F*invV-R)", s'il vous plaît ?

attention c'est plutot det(F-r*V) que det(V-r*F)
car (F*invV-R)*V=F-r*V donc si det(V) est non nul vous avez

det(F*invV-R)=det(F-r*V) a facteur non nul pres, vous aurez donc les memes racines en r


vous etes en mode maple? car sinon les tableaux comme Vbase commencent a 0 ca doit etre pour cela que pour moi V avait un det nul.

[Legion]

Non, et meme si la factorisation aboutit, ca donnerait surement des resultats beaucoup trop gros donc inutilisables au pire et tres inefficaces au mieux. Je pense que l'approche vouloir une solution analytique n'est pas la bonne (un peu comme vouloir une formule pour les equations du 4eme degre). Il serait sans doute beaucoup plus efficace de remplacer le script actuel par une fonction prenant en argument les valeurs des parametres et renvoyant le resultat du solve(...) dans ce cas.
Donc remplacer le mot clef purge par f puis mettre une accolade { apres la parenthese fermante a la place du point virgule et terminer par return res; }

Je suis en mode maple.
Je comprends que je suis sur la mauvaise voie. Je vais essayer de voir comment faire le lien entre xcas et mon programme en R.
Merci beaucoup à vous tous pour votre aide.