Discussion :

[keaton84]

Hello tout le monde,


Comme le titre de mon poste l'indique, j'aimerais calculer le résultat d'un dév limité d'une expo matricielle, le truc c'est que je ne veux pas utiliser de boucle...

la formule d'un dév limité d'une expo est assez simple :
Pour P une matrice,

Exp(P) = Somme (P^i/ factoriel (i)), pour i --> à l'infini

Je sais aussi que P converge vers une matrice nulle (on dit que P est nilpotente), du coup mon i n'ira (heureusement) jamais vers l'infini, à condition bien sûr que je lui trouve une jolie condition d'arrêt dans mon code

Il est assez facile de faire ce calcul en utilisant une boucle, mais si je veux profiter de la structure matricielle de R, y a surement moyen de le faire via un vecteur c(1:infini) , dont chaque élément élèvera P à la puissance i .

Je bloque sur deux choses:

> avec quel symbole/fonction dans R, je peux faire : sum(P^i) sachant que i c'est les éléments de c(1:infini)
> comment introduire une condition d'arrêt au niveau d'un entier q tel que P^q < matrice Epsilon

Merci par avance de votre aide (Ma demande n'est peut être pas super claire mais y a qu'avec vos remarques que je le saurai )

Réda

[mamounMob]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
lapply(1:10), function(x) {
                    Mat^x
}

'lapply' te renvoie une liste qui contiendra les puissances de ta matrice, tu pourra faire la somme après et poser ta condition d'arrêt.

[keaton84]

Hello Mamoun, sympa merci

Finalement la fonction Expo Matricielle existe déjà, c'est expm() (ça s'invente pas lol), elle approche la série infinie que je décris plus haut via un algo d'approximation. voici ce que je trouve :

The exponential of a matrix is defined as the infinite Taylor series expm(A) = I + A + A^2/2! + A^3/3! + ... (although this is definitely not the way to compute it). The method for the dgeMatrix class uses Ward's diagonal Pade' approximation with three step preconditioning.

A++

[keaton84]

Hello,

J'ai appliqué la solution de Mamoun mais j'ai l'impression que m^x me calcule chaque élément de P à la puissance x alors que moi c'est la puissance matricielle dans le sens produit matricielle qui m'intéresse, càd:

P^1=P
P^2=P%*%P
P^3=P%*%P%*%P
....

Une idée sur la syntaxe ?


Merci

[mamounMob]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
require(expm)
mat=matric(1:4, 2)
mat%^%2

Voilà faut juste ajouter les '%'

[keaton84]

Merci mec, ça marche nickel !!

En revanche, je bloque sur un truc maintenant, j'ai le code suivant :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
P<-(Matrix_Transition-diag(12))
N<-10
Serie<-lapply(1:N,function(u) P^u)

Qui me donne:

Serie[1] = P
Serie[2] = P^P
Serie[3] = P^P^P
...

Je veux sommer Serie[1]+Serie[2]+....+Serie[N]=P_Globale

P_Globale étant une matrice de même taille que P mais dont les éléments représentent la somme des N premières matrices...

P_Globale(1,1) = Serie[1](1,1)+Serie[2](1,1)+....+Serie[N](1,1)
P_Globale(2,1) = Serie[1](2,1)+Serie[2](2,1)+....+Serie[N](2,1)
...

Sauriez-vous comment faire ? Sachant que Serie[i] n'est même pas considérée comme numérique

Merci d'avance.