Discussion :

[saad.pizzicato]

Bonjour tout le monde,

Ma question concerne la résolution des systèmes d'équations algébriques du type f(x)=0 ou f={fi} un vecteur de n fonctions, x=(xi) le vecteur des inconnues...chacune des fi, étant une fonction non linéaire des xi...

Ma question est si MATLAB est déjà doté des algorithmes de résolution de ce genre d'équations (Newton Raphson, deepest descent...) ou que je dois moi même programmer mon algorithme de résolution???


Merci à l'avance...

[shaiHulud]

Bonsoir,

Pour les descentes de gradients, doc fminunc et fmincon. (Pour trouver un zero, on peut minimiser f(x)^2).

[saad.pizzicato]

Bonsoir,

Merci pour vos indications;

J'ai lu sur la doc de MATLAB que la commande fsolve:

fsolve finds a root (zero) of a system of nonlinear equations.

Sauf que j'ai trouvé nulle part les algorithme utilisés par le solveur MATLAB.

A part cela, et pour tester cette commande: comment je peux définir mes fonctions pour l'utiliser?
Doit je les définir dans un fichier à part? ou je peux la définir directement...Ma fonction est la suivante:

Code :

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

1
2
3
r=2:n;
a(1)=0,05;
f(r)=2*sum(a(i)*k(i,r),i)+3*sum(sum(sum(a(i)*a(j)*a(k)*b(i,j,k,r),k),j),i)-2*[{sum(sum(a(i)*a(j)*k(i,j),j),i)+sum(sum(sum(sum(a(i)*a(j)*a(k)*a(j)*a(l)*b(i,j,k,l),l),k),j),i)}/(sum(sum(a(i)*a(j)*m(i,j),j),i))]*sum(a(i)*m(i,r),i) = 0

Ou je résous les a(i) (i=1,n) pour les f(r) = 0, avec r=2:n et le a(1) est une valeur que j'injecte ) chaque fois, c'est à dire résoudre un système de n-1 équation à n variable ou j'ai le choix de fixer cette variable a(1)...

Si c'est pas claire mon système, le voilà sur l'image attachée, ça va mieux apparaitre (N.B: les sommation sont écrite par convention de sommation d'Einstein).

Merci à l'avance...

[le fab]

salut

il faut écrire une et une seule fonction qui décrit ton système, de la forme

Code :

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

function f = myfun(a)

ou f est un vecteur de taille n-1 et a un vecteur de taille n

ensuite dans l'appel de fsolve tu peux utiliser @ pour fixer a(1) de la sorte :

Code :

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

fsolve(@(x) myfun([valA1 x],....)(

ou x représente tes inconnues a(2:end)

Fabien

[saad.pizzicato]

Bonjour,
Je bloque toujours sur ce qui suit:
Juste pour essai, j'ai créé un fichier nommé SystAlg qui comporte la définition suivante:

Code :

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

1
2
3
4
5
6
7
8
9
 
function f = SystAlg (a)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
 
for r=2:25;
f(a,r)= 2*sum(a(i)*K(i,r),1);
end
end

Puis dans le code initial j'ai tapé:

Code :

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

1
2
3
 
a(1)=0.05;
a=fsolve(@a(1)SystAlg,a(2:25));

Sauf que je reçois l'erreur suivante:

??? Error: File: Plaque.m Line: 166 Column: 10
Unbalanced or unexpected parenthesis or bracket.

??? Error: File: Plaque.m Line: 166 Column: 12
Unbalanced or unexpected parenthesis or bracket.

Merci à l'avance...

[saad.pizzicato]

Une autre question:

Comment je peux écrire sur Matlab:

Sigma(i,Sigma(j, Sigma (k, ai.aj.ak.bijkr))) = Triple sommation selon i, j et k de ai.aj.ak.b(i,j,k,l).


Merci à l'avance...

[tachmou]

Bonjour,

Le code que vous avez posté n'est pas complet... je n'arrive donc pas à comprendre ce que vous cherchez à faire exactement.

Cependant, voici un exemple qui se rapproche le plus de ce que vous cherchez à faire. Enregistrez ce code dans un fichier que vous nommerez "test.m" et exécutez-le. J'espère que cela va vous aider un peu :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
function test
 
    a(1) = 0.05;                 %abscisse de départ pour la recherche
 
    a_optimized = fsolve(@f, a)  %abscisse trouvée par la minimisation
 
    X = [-1:0.01:1];
    Y = f(X);
 
    plot(X,Y)
    hold on
    grid on
    xlim([-1 1]);
    ylim([-1 1]);
    plot(a_optimized, 0, 'or');  % Affichage du point trouvé avec un cercle (o) rouge (r)
 
    function out = f(x)
        out = x.*x;
    end
end

La fonction fsolve a besoin de l'adresse de la fonction a minimiser : grâce à cette adresse, fsolve peut l'utiliser comme bon lui semble pour effectuer sa minimisation. C'est ce qu'on appelle un handle dans matlab. On l'obtient via le symbole "@".

Pour effectuer votre somme, vous pouvez procéder comme suit :

Code :

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

1
2
3
4
clear M %efface ce que contient M
M(2,2,3) = 0 %crée une matrice 3x3x3 remplie de 0
M(2,1,3) = 2.24 % affecte la valeur 2,24 à la coordonnée 2,1,3
somme_totale = sum(sum(sum(M)))

Il vous suffit donc de créer une matrice M contenant le produit des éléments qui vous intéresse. Ici, M(i,j,k) = a(i)*a(j)*a(k)*b(i,j,k,l) pour tous les i,j,k. Cela peut être la sortie d'une fonction de l pour traiter le paramètre l. Il faut aussi que b(i,j,k,l) soit définie comme une fonction, mais de i,j,k,l cette fois-ci.

N'hésitez pas à modifier l'exemple "test.m" pour bien comprendre comment tout ça fonctionne. N'hésitez pas à poser d'autres questions.

Cordialement,

[saad.pizzicato]

Bonjour,
en fait j'ai définis le système de fonctions moyennant la méthode matricielle que vous avez proposé, pour l'instant c'est ça ce que ça donne:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
function f = myfun(a)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

Clear a
a(1)=0.05;  %Je fixe ma valeur de a(1)%
Term1(25,25)=0;
for r=1:25
    for i=1:25
        Term1(i,r)=a(i)*K(i,r);
    end
end

TERME1=sum(Term1,1);

Term2(25,25,25,25)=0;
for r=1:25
    for i=1:25
        for j=1:25
            for k=1:25
                Term2(i,j,k,r)=a(i)*a(j)*a(k)*B(i,j,k,r);
            end
        end
    end
end

TERME2=sum(Term2,4);

AK(25,25)=0;
for i=1:25
    for j=1:25
        AK(i,j)=a(i)*a(j)*K(i,j);
    end
end

SOMAK=sum(sum(AK));

AB(25,25,25,25)=0;
for i=1:25
    for j=1:25
        for k=1:25
            for l=1:25
        AB(i,j,k,l)=a(i)*a(j)*a(k)*a(l)*B(i,j,k,l);
            end
        end
    end
end
    
SOMAB=sum(sum(sum(sum(AB))));

AM(25,25)=0;
for i=1:25
    for j=1:25
        AM(i,j)=a(i)*a(j)*M(i,j);
    end
end

SOMAM=sum(sum(AM));

Term3(25,25)=0;
for r=1:25
    for i=1:25
        Term3(i,r)=a(i)*M(i,r);
    end
end

TERME3=sum(term3,1);

for r=2:25
    f(r)=(2*TERME1(r))+(3*TERME2(r))-(2*((SOMAK+SOMAB)/SOMAM)*TERME3(r));
end   %Ma fonction à minimiser%
    
end

Votre avis tout d'abord sur la saisie de la fonction...
Et

Il y a un éclaircissement à faire:
C'est que étant donné a=a(i) le vecteur des inconnues...le a(1)=0,05 n'est pas une valeur de départ...Ce que je veux c'est FIXER a(1) à cette valeur et retrouver, en fonction de la valeur FIXEE pour a(1) les autres a(i)....
Parce que comme vous avez bien pu remarquer: il y a n-1 équations pour n inconnues...donc je fixe a(1) pour n'en avoir que n-1 et ça devient faisable...
Alors comment traduire cette donnée??

Ensuite j'ai évoqué ma fonction dans le programme initial comme tel:

Code :

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

1
2
3
4
 
 
 
a= fsolve(@myfun)

Est ce la bonne démarche???
Merci à l'avance...

[le fab]

salut

tu ne t'y prend pas bien

déjà matlab est capable de résoudre un système de n-1 équations de n inconnues ... certes ce ne sera pas une solution unique mais il te proposera une solution. exemples :

Code :

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

1
2
3
fsolve(@(x) x(1)+x(2),[0 0])
fsolve(@(x) x(1)+x(2),[1 0])
fsolve(@(x) x(1)+x(2),[-3 12])

ensuite fsolve prend deux arguments (au moins)
la fonction et le point de départ de l'optimisation

Code :

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

fsolve(myfun,a0)

surtout, le clear a en début de fonction est à éviter absolument, surtout que là, a est le paramètre d'entrée de ta fonction

enfin si tu veux fixer a(1) il vaut mieux le faire en dehors de la fonction et appeler dans fsolve une fonction à n-1 inconnues, via l'operateur @ par exemple, comme suit :

Code :

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

1
2
3
4
5
 
a1 = 0.05;
a2_end_0 = zeros(n,1); % un exemple de point de départ, pas forcement pertinent
a2_end = fsolve(@(x) myfun([a1;x]),a2_end_0 );
a = [a1;a2_end ];

en prenant soin de virer la ligne a1 = 0.05; au début de ta fonction

Fabien

[tachmou]

Bonjour,

Désolé de répondre aussi tard, mais j'ai été très pris ces derniers temps.

Ne prenez pas ce que je vais dire comme un affront, mais comme un conseil à suivre qui va vous permettre de répondre vous même à vos questions.

Tout d'abord, j'ai vu dans votre code que vous avez compris beaucoup de choses. C'est excellent.

Par ailleurs, j'approuve les commentaires de "le fab". Je vous conseille de méditer dessus.

Maintenant, prenons un peu de recul.

Vous vous attachez à un problème un peu complexe (vous avez bien vu que le nombre de lignes de codes a vite augmenté). Je vous donne un premier conseil d'ordre général :
1) J'ai cru comprendre que vous étiez enseignant-chercheur. C'est là que vous risquez de prendre mon conseil pour un affront, mais il ne faut pas. La méthode de base utilisée pour faire une recherche en général est d'écrire sur une feuille en haut ce que l'on a et en bas ce que l'on cherche, et d'essayer de raccorder les deux bouts via des étapes intermédiaires. Je vous conseille d'appliquer ce principe pour les étapes générales de la résolution, et pour chaque étape, ré appliquer le même principe pour détailler la démarche. Arrivé à 3 ou 4 niveau de détails, on commence à avoir une bonne vue.
2) Une fois que vous avez appliqué le 1), vous voyez de quoi vous avez besoin. Je vous conseille donc de programmer des petites fonctions qui feront les calculs "de base". Puis vous programmer les fonctions du niveau "supérieur" en utilisant les fonctions du niveau "de base". Et ainsi de suite en allant jusqu'au programme principal.

Je vous conseille sincèrement d'essayer ce que je vous dit, en essayant d'avoir un bon recul tout le long de la résolution. Cela a fait ses preuves chez plusieurs personnes : ça structure la pensée et ça aide fortement à la résolution. De plus, cet exercice, en l'appliquant correctement sur matlab, vous fera faire de très bons et rapides progrès pour maîtriser les calculs "élémentaires". Enfin, la programmation des étapes de base consiste simplement en la programmation de fonctions possédant des entrées (ce que vous avez) et des sorties (ce que vous cherchez).

Si je vous dit ça, c'est que j'ai bien constaté que vous étiez capable de le faire. Je veux bien vous aider pour chacune des étapes de base, et on les valide ensemble. Il sera ainsi d'autant plus simple de vous corriger, car là, ça reste un peu complexe (et surtout long car là, je manque de temps).

Enfin, le fait de procéder ainsi permet aussi de valider les étapes les unes après les autres en validant les sous-fonctions, ce qui est plus simple que de valider un programme entier d'un seul coup.

Après, s'il y a trop de "sous-fonctions", il y a toujours la possibilité de simplifier par la suite si vous voulez.

Cordialement,

[saad.pizzicato]

Bonsoir tout le monde,

Je suis entrain d'essayer tout d'abord la question évoquée par Le FAB:

résoudre un système de n-1 équations de n inconnues ... certes ce ne sera pas une solution unique mais il te proposera une solution. exemples :

Après dès résolution, je passera à mon problème réel.

Ma fonction est maintenant:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
function f = myfun(a)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
 
M=load('M.mat');
K=load('K.mat');
B=load('B.mat');
 
%J'ai ajouté ces lignes puisqu'au début je recevai un message d'erreur comme quoi il ne reconnais pas mon K, B, M du programme de base, ces matrices sont stockée dans des fichiers%
 

 
Term1(25,25)=0;
for r=1:25
    for i=1:25
        Term1(i,r)=a(i)*K(i,r);
    end
end
 
TERME1=sum(Term1,1);
 
Term2(25,25,25,25)=0;
for r=1:25
    for i=1:25
        for j=1:25
            for k=1:25
                Term2(i,j,k,r)=a(i)*a(j)*a(k)*B(i,j,k,r);
            end
        end
    end
end
 
TERME2=sum(Term2,4);
 
AK(25,25)=0;
for i=1:25
    for j=1:25
        AK(i,j)=a(i)*a(j)*K(i,j);
    end
end
 
SOMAK=sum(sum(AK));
 
AB(25,25,25,25)=0;
for i=1:25
    for j=1:25
        for k=1:25
            for l=1:25
        AB(i,j,k,l)=a(i)*a(j)*a(k)*a(l)*B(i,j,k,l);
            end
        end
    end
end
    
SOMAB=sum(sum(sum(sum(AB))));
 
AM(25,25)=0;
for i=1:25
    for j=1:25
        AM(i,j)=a(i)*a(j)*M(i,j);
    end
end
 
SOMAM=sum(sum(AM));
 
Term3(25,25)=0;
for r=1:25
    for i=1:25
        Term3(i,r)=a(i)*M(i,r);
    end
end
 
TERME3=sum(term3,1);
 
for r=2:25
    f(r)=(2*TERME1(r))+(3*TERME2(r))-(2*((SOMAK+SOMAB)/SOMAM)*TERME3(r));
end
    
end

J'ai évoqué ma fonction dans le programme de base par:

Code :

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

1
2
3
4
5
6
7
8
9
a0(1)=0.05;
for i=2:25;
    a0(i)=0;
end
 
%Point de départ pour la résolution%

a=fsolve(@myfun,a0,optimset('Display','iter'));

Je reçois l'erreur suivante:

??? Undefined function or method 'mtimes' for input
arguments of type 'struct'.

Error in ==> myfun at 15
        Term1(i,r)=a(i)*K(i,r);

Error in ==> fsolve at 248
            fuser = feval(funfcn{3},x,varargin{:});

Error in ==> Plaque at 178
a=fsolve(@myfun,a0,optimset('Display','iter'));

Caused by:
    Failure in initial user-supplied objective function
    evaluation. FSOLVE cannot continue.

Merci à l'avance...

PS: Pour M.Tachmou, un petit éclaircissement : je ne suis pas enseignant chercheur mais plutôt ingénieur chercheur au tout début de ma thèse... je prendrai jamais ce que vous dite comme affront, au contraire, tout conseil est le bienvenu...Finalement on apprend et on s'améliore avec le temps. Mais je tiens quand même à préciser que je suis pas vraiment loin de ce que vous dites puisque déjà j'ai parcouru une bonne distance dans l'élaboration de la partie théorie, pour l'application numérique je me plante puisque j'avoue que l'informatique et la programmation c'était une terre inconnue pour moi (Parfois dans la formulation, parfois dans la synthaxe)...J'apprends par vos conseils...Merci encore une fois.

[le fab]

salut

que te retourne

Code :

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

whos K

il semblerait qu'il s'agisse d'une structure / tableau de structure, donc il doit te manquer quelque chose comme :

Code :

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

Term1(i,r)=a(i)*K(i,r).champ;

sinon

Code :

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

1
2
3
4
a0(1)=0.05;
for i=2:25;
    a0(i)=0;
end

s'écrit plus simplement

Code :

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

1
2
a0 = zeros(1,25);
a0(1)=0.05;

ou

Code :

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

a0 = [0.05 zeros(1,24)];

[saad.pizzicato]

Bonjour,

J'ai pu, soit-disant résoudre le problème en mettant tout simplement load('M'); au lieu de M=load('M'); dans la définition de ma fonction...

Apparemment, le programme commence la résolution de mon système, ce qui est déjà pas mal pour moi Mais:

les problème que je rencontre maintenant sont les suivants:

Cas n°1:

En essayant, de résoudre le problème de n inconnues à n-1 équation j'ai ceci:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
>> a=fsolve(@myfun,a0,optimset('Display','iter'));
 
                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0         26    3.03726e-025                     6.08e-024               1
 
Equation solved at initial point.
 
fsolve completed because the vector of function values at the initial point
is near zero as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.
 
<stopping criteria details>

J'ai essayé d’agrandir la tolérance d'une manière un peu arbitraire pour voir un peu l'effet que ça donne:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
>>  a=fsolve(@myfun,a0,optimset('Display','iter','TolFun',1e-40));
 
                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0         26    3.03726e-025                     6.08e-024               1
     1         27    3.03726e-025              0      6.08e-024               1
 
No solution found.
 
fsolve stopped because the problem appears to be locally singular.
 
<stopping criteria details>

Cas n°2:

Qui est mon problème de résoudre ce système pour a(1)=0.05 et calculer les a(2:25):

Code :

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

1
2
3
4
5
%Dans la définition de ma fonction j'ai ajouté l'argument x%
function f = myfun(x,a)
 
%Par la suite dans la définition de la fonction%
x=a(2:25);

Dans le code principal:

Code :

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

1
2
3
a0(2:25)=1/1000;
>> a(1)=0.05;
>> a=fsolve(@(x) @myfun,a0,optimset('Display','iter'));

Et l'erreur

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using
Levenberg-Marquardt algorithm instead. 
> In fsolve at 314
??? Undefined function or method 'isfinite' for input arguments of type 'function_handle'.

Error in ==> levenbergMarquardt at 14
if any(~isfinite(costFun))

C'est correct ce que j'ai fais dans les deux cas? Sinon les problèmes susmentionnés viennent d’où à votre avis?


Mes salutations...

[le fab]

Bonjour,

J'ai pu, soit-disant résoudre le problème en mettant tout simplement load('M'); au lieu de M=load('M'); dans la définition de ma fonction...

alors la première syntaxe charges les variables de M dans le workspace, la seconde dans une structure M, d'ou ton erreur

Cas n°2:

Qui est mon problème de résoudre ce système pour a(1)=0.05 et calculer les a(2:25):

Code :

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

1
2
3
4
5
%Dans la définition de ma fonction j'ai ajouté l'argument x%
function f = myfun(x,a)
 
%Par la suite dans la définition de la fonction%
x=a(2:25);

ca me parait foireux !

Dans le code principal:

Code :

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

1
2
3
a0(2:25)=1/1000;
>> a(1)=0.05;
>> a=fsolve(@(x) @myfun,a0,optimset('Display','iter'));

plutot quelque chose comme

Code :

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

a=fsolve(@(x) myfun(x,a),a0,optimset('Display','iter'));

regarde bien la doc de fsolve, j'ai l'impression que tu te mélange un peu sur la syntaxe
c'est d'ailleur la principale difficulté en optimisation sous matlab : écrire son problème au format attendu par matlab

[saad.pizzicato]

Merci beaucoup, ça fonctionne maintenant...

Sinon, je remarque que la vitesse de convergence est très très lente...Y a t-il des méthode pour accélérer la convergence pour fsolve? des paramètres à modifier etc....

[le fab]

salut

c'est quoi "très très lent" ?

sinon pour accelerer tu peux :
- définir un point de départ x0 plus pertinent (plus proche de la solution) si c'est possible
- changer d'algorithme
- utiliser (si possible) un jacobian (jacobian, jacomult, jacopattern). mais bon, c'est pas toujours possible / facile de le faire

Fabien

[saad.pizzicato]

Merci beaucoup,
En fait j'ai changé l'algorithme en adoptant celui 'levenberg-marquardt' et ça me donne exactement ce que je veux...

Merci encore une fois...

[le fab]

alors