Bonjours a tous.
Voici le système d'équation correspondant a la version a
une dimension du problème physique.
A0:matrix(
[ w1, w2, w1, 0, 0, 0, 0, 0, 0, 0, 0, p1],
[ 0, w1, w2, w1, 0, 0, 0, 0, 0, 0, 0, p2],
[ 0, 0, w1, w2, w1, 0, 0, 0, 0, 0, 0, p3],
[ 0, 0, 0, w1, w2, w1, 0, 0, 0, 0, 0, p4],
[ 0, 0, 0, 0, w1, w2, w1, 0, 0, 0, 0, p5],
[ 0, 0, 0, 0, 0, w1, w2, w1, 0, 0, 0, p6],
[ 0, 0, 0, 0, 0, 0, w1, w2, w1, 0, 0, p7],
[ 0, 0, 0, 0, 0, 0, 0, w1, w2, w1, 0, p8],
[ 0, 0, 0, 0, 0, 0, 0, 0, w1, w2, w1, p9],
[ w3, w4, w5, w4, w3, 0, 0, 0, 0, 0, 0, pA],
[ 0, w3, w4, w5, w4, w3, 0, 0, 0, 0, 0, pB],
[ 0, 0, w3, w4, w5, w4, w3, 0, 0, 0, 0, pC],
[ 0, 0, 0, w3, w4, w5, w4, w3, 0, 0, 0, pD],
[ 0, 0, 0, 0, w3, w4, w5, w4, w3, 0, 0, pE],
[ 0, 0, 0, 0, 0, w3, w4, w5, w4, w3, 0, pF],
[ 0, 0, 0, 0, 0, 0, w3, w4, w5, w4, w3, pG]
);
L'un des problèmes qui empêche l'utilisation des fonctions
Maxima est que le système a plus d'équations que d'inconnues
Il y a donc des équations redondantes.
J'ai essayé de résoudre ce petit système avec mon code mais même
dans ce cas la capacité mémoire utilisé par maxima grossie de façon
inconsidéré. Le phénomène se produit au moment de l'appel de la fonction
"fullratsimp()" sur une seul expression. J'ai essaye avec "ratsimp()" mais
le résultat est identique. La mémoire dépasse les 200Mo alors qu'elle est
à moins de 40Mo avant l'appel de cette fonction.
Voici le code que j'utilise actuellement. La fonction "complexity()"
mériterait d'être optimisé.
MODE_DECLARE(FUNCTION(COMPLEXITY),FIXNUM);
/*
* Calcule le nombre d'atome dans la forme
* EXPANDed de l'expression.
*
*/
COMPLEXITY(EXPR):=BLOCK(
[I,M:0,N],
MODE_DECLARE([I,M,N],FIXNUM),
IF EXPR=0 THEN RETURN(0),
IF ATOM(EXPR) THEN RETURN(1),
IF OP(EXPR)="-" THEN EXPR:-EXPR,
IF ATOM(EXPR) THEN RETURN(1),
N:LENGTH(EXPR),
IF OP(EXPR)="+" THEN (
FOR I:1 THRU N DO (
M:M+COMPLEXITY(PART(EXPR,I))
)
) ELSE (
EXPR:EXPAND(EXPR),
IF OP(EXPR)="*" THEN RETURN(N),
IF OP(EXPR)="^" THEN RETURN(N-1),
M:COMPLEXITY(EXPR)
),
EXPR:FALSE,
M
);
EVAL_WHEN(TRANSLATE,DECLARE_TRANSLATED(COMPLEXITY));
SOLVEXI(A0,NCOL,LASTROW,PIVOTS):=BLOCK(
[NROW,I,J,K,L,K1,J1,CI,NCI,CJ,NCJ,KI,KJ,MINDEX,M,MVAL,N,S,PGCD,PIVOT,LNT,LSNT:[],LP,LPGCD:[],ML
PGCD:[]],
MODE_DECLARE([NROW,NCOL,LASTROW,I,J,K,L,K1,J1,NCI,NCJ,MINDEX,M,MVAL,PIVOT,N],FIXNUM),
NROW:LENGTH(A0),
N:0,
FOR I:1 THRU NROW DO (
S:[],
M:0,
LPGCD:CONS(1,LPGCD),
MLPGCD:CONS(1,MLPGCD),
FOR J:NCOL+1 STEP -1 THRU 1 DO (
IF A0[I,J]=0 THEN (
S:CONS(0,S)
) ELSE (
IF ABS(A0[I,J])=1 THEN (
S:CONS(1,S),
M:M+1
) ELSE (
N:COMPLEXITY(A0[I,J]),
S:CONS(N,S),
M:M+N
)
)
),
LSNT:ENDCONS(M,LSNT),
IF I=1 THEN (
LNT:MATRIX(S)
) ELSE (
LNT:ADDROW(LNT,S)
)
),
FOR I:1 THRU LASTROW DO (
MVAL:2147483648,
PIVOT:-1,
/* RECHERCHE DU MEILLEUR PIVOT
*
* POUR CHAQUE COLONNE SAUF LA DERNIERE
* RECHERCHER LE MEILLEUR PIVOT DE CETTE COLONNE
*/
FOR J:1 THRU NCOL DO (
IF PIVOTS[J] THEN (
N:0,
S:FALSE,
LP:[],
/* CONTRUIT LA LISTE DES INDICES CORRESPONDANT
* A DES CELLULE NON NUL DE LA MATRICE DANS LA
* COLONNE J, ET DETERMINE SI L'UNE D'ENTRE ELLE
* AU MOINS PEUT ETRE PRISE COMME PIVOT.
*/
FOR K:1 THRU NROW DO (
IF LNT[K,J]#0 THEN (
N:N+1,
LP:CONS(K,LP),
IF PIVOTS[K] THEN S:TRUE
)
),
/* SI LE NOMBRE DE CELLULE NON VIDE DE LA MATRICE
* EST > 0 ET QUE AU MOINS UNE D'ENTRE ELLE PEUT
* ETRE PRISE COMME PIVOT ALORS POUR CHAQUE PIVOT
* POTENTIEL DE CETTE COLONNE ESTIMER AU MIEUX
* L'AUGMENTATION DE TAILLE DU SYSTEME RESULTANT.
*/
IF N>0 AND S THEN (
/* POUR CHAQUE PIVOT POTENTIEL PRECEDEMENT
* DETERMINER.
*/
FOR K:1 THRU N DO (
/* RECUPERER L'INDEX DE CETTE CELLULE */
K1:LP[K],
/* SI CETTE CELLULE PEUT ETRE PRISE COMME PIVOT
* ALORS CALCULER L'AUGMENTATION DE TAILLE DU
* SYSTEME RESULTANT
*/
IF PIVOTS[K1] THEN (
M:0,
CI:A0[K1,J],
FOR L:1 THRU N DO (
J1:LP[L],
IF L#K THEN (
PGCD:GCD(CI,A0[J1,J]),
IF ABS(PGCD)#1 THEN (
LPGCD[J1]:PGCD,
NCI:COMPLEXITY(RATSIMP(CI/PGCD)),
NCJ:COMPLEXITY(RATSIMP(A0[J1,J]/PGCD))
) ELSE (
NCI:LNT[K1,J],
NCJ:LNT[J1,J]
),
M:M+NCI*LSNT[J1]+NCJ*LSNT[K1],
M:M-NCI*LNT[J1,J]-NCJ*LNT[K1,J],
IF M>MVAL THEN RETURN(0)
)
),
IF N=1 THEN M:LSNT[K1],
IF M9 AND J>12 AND K=12 THEN (
A0[J,K]:FULLRATSIMP(EXPAND(KJ*A0[PIVOT,K]-KI*A0[J,K]))
) ELSE (
A0[J,K]:FULLRATSIMP(KJ*A0[PIVOT,K]-KI*A0[J,K])
),
LNT[J,K]:COMPLEXITY(A0[J,K]),
IF PGCD=0 THEN (
PGCD:A0[J,K]
) ELSE (
IF ABS(PGCD)#1 THEN (
PGCD:GCD(PGCD,A0[J,K])
)
)
),
FOR K:PIVOT+1 THRU NCOL+1 DO (
/* a moins que ce ne soit ici */
IF I>9 AND J>12 AND K=12 THEN (
A0[J,K]:FULLRATSIMP(EXPAND(KJ*A0[PIVOT,K]-KI*A0[J,K]))
) ELSE (
A0[J,K]:FULLRATSIMP(KJ*A0[PIVOT,K]-KI*A0[J,K])
),
LNT[J,K]:COMPLEXITY(A0[J,K]),
IF PGCD=0 THEN (
PGCD:A0[J,K]
) ELSE (
IF ABS(PGCD)#1 THEN (
PGCD:GCD(PGCD,A0[J,K])
)
)
),
KI:0,
KJ:0,
IF ABS(PGCD)#1 THEN (
IF PGCD#A0[J,NCOL+1] THEN (
DISPLAY(PGCD),
M:0,
FOR K:1 THRU NCOL DO (
IF I>9 AND J>12 THEN DISPLAY(K),
A0[J,K]:FULLRATSIMP(A0[J,K]/PGCD),
LNT[J,K]:COMPLEXITY(A0[J,K]),
M:M+LNT[J,K]
),
LSNT[J]:M
) ELSE (
LSNT[J]:COMPLEXITY(A0[J,NCOL+1])
)
),
PGCD:0
)
)
),
DISPLAY(I),
CI:0
),
[A0,PIVOTS]
);
dynamalloc:true;
facexpand:false;
AVAILABLE_PIVOTS_LIST:[];
for i:1 thru 16 do (
AVAILABLE_PIVOTS_LIST:CONS(TRUE,AVAILABLE_PIVOTS_LIST)
);
result:SOLVEXI(A0,11,11,AVAILABLE_PIVOTS_LIST);
Laurent