résolution d'un systeme non linéaire

**hichamx** · 15/05/2011, 00h40

Salut à vous tous
bon je suis sur la réalisation d'un mini projet d'un programme de résolution d'un système non linéaire à 3 équations et 3 variables par la méthode de Newton - Raphson
le problème c'est que si j'initialise les x0, y0 et z0 par des nombres qui différent de la solution par 0.01 par exemple alors la solution diverge
en plus chaque fois si je veux afficher le nombre d’itérations (compt) il m'affiche toujours 1
voici le code

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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
 
float f (float x,float y,float z)
{ return (x*x+y*y+z*z-9);}
float g (float x,float y,float z)
{ return (x*x+y*y-3*x);}
float h (float x,float y,float z)
{ return (y+z-2);}
float df_1(float x,float y,float z)
{ return (2*x);}
float df_2(float x,float y,float z)
{ return (2*y);}
float df_3(float x,float y,float z)
{ return (2*z);}
float dg_1(float x,float y,float z)
{ return (2*x-3);}
float dg_2(float x,float y,float z)
{ return (2*y);}
float dg_3(float x,float y,float z)
{ return (0);}
float dh_1(float x,float y,float z)
{ return (0);}
float dh_2(float x,float y,float z)
{ return (1);}
float dh_3(float x,float y,float z)
{ return (1);}
 
 
int main()
{
 
 
float tab[3][4];
float tab2[3][4];
int i,j,compt;
float x,y,z,x0,y0,z0;
float  epsilon = pow(10,-1000000); 
x0=2.645;
y0=0.9685;
z0=1.0314;
 
 
 
 
    tab[0][0]= df_1( x0,y0,z0);
    tab[0][1]= df_2( x0,y0,z0);
    tab[0][2]= df_3( x0,y0,z0);
    tab[0][3]= (df_1( x0,y0,z0))*x0+(df_2( x0,y0,z0))*y0+(df_3( x0,y0,z0))*z0-f ( x0,y0,z0);
    tab[1][0]= dg_1( x0,y0,z0);
    tab[1][1]= dg_2( x0,y0,z0);
    tab[1][2]= dg_3( x0,y0,z0);
    tab[1][3]= (dg_1( x0,y0,z0))*x0+(dg_2( x0,y0,z0))*y0+(dg_3( x0,y0,z0))*z0-g ( x0,y0,z0);
    tab[2][0]= dh_1( x0,y0,z0);
    tab[2][1]= dh_2( x0,y0,z0);
    tab[2][2]= dh_3( x0,y0,z0);
    tab[2][3]= (dh_1( x0,y0,z0))*x0+(dh_2( x0,y0,z0))*y0+(dh_3( x0,y0,z0))*z0-h ( x0,y0,z0);
 
 
 
	for (i=0;i<4;i++)		
	{
		tab2[0][i]=tab[0][i];
	}
 
 
 
	for (i=0;i<4;i++)		
	{
		tab2[1][i]=(tab[1][i]*tab[0][0])-(tab[0][i]*tab[1][0]); 
 
	}
 
	for (i=0;i<4;i++)
	{
		tab2[2][i]=(tab[2][i]*tab[0][0])-(tab[0][i]*tab[2][0]);
 
	}
 
 
 
for (i=0;i<3;i++)
{
	printf("\n\n");
 
	for(j=0;j<4;j++)
	{
		printf("\t %f",tab2[i][j]);
	}
}
 
		for (i=0;i<3;i++)	
		{
			for (j=0;j<4;j++)
			{
				tab[i][j]=tab2[i][j];
			}
 
		}
 
 
 
	for (i=0;i<2;i++)
	{
		for (j=0;j<4;j++)
		{
			tab2[i][j]=tab[i][j];
		}
	}
 
 
	tab2[2][1]=(tab[1][1]*tab[2][1])-(tab[1][1]*tab[2][1]);
	tab2[2][2]=(tab[1][2]*tab[2][1])-(tab[1][1]*tab[2][2]);
	tab2[2][3]=(tab[1][3]*tab[2][1])-(tab[1][1]*tab[2][3]);
 
 
 
 
 
for (i=0;i<3;i++)
		{
			printf("\n\n\n");
 
			for (j=0;j<4;j++)
			{
				printf("\t %f ",tab2[i][j]);
			}
 
		}
 
		z=tab2[2][3] / tab2[2][2];
		y=((tab2[1][3]-(tab2[1][2]*z)))/tab2[1][1];
		x=((tab2[0][3]-(tab2[0][2]*z)-(tab2[0][1]*y)))/tab2[0][0];
 
compt=1;
 
while( abs(x-x0)>epsilon  &&  abs(y-y0)>epsilon  &&  abs(z-z0)>epsilon )
 
{		
        x0=x;
        y0=y;
        z0=z;
 
        z=tab2[2][3] / tab2[2][2];
		y=((tab2[1][3]-(tab2[1][2]*z)))/tab2[1][1];
		x=((tab2[0][3]-(tab2[0][2]*z)-(tab2[0][1]*y)))/tab2[0][0];
 
        compt=compt+1;
 
 
 
 
}
 
printf("\n\n x =%f",x);
printf("\n   y =%f",y);
printf("\n	 z =%f \n\n",z);
printf("le nombre d'iterations est %d\n",compt);
system("PAUSE");
}

**Sve@r** · 15/05/2011, 20h07

Envoyé par hichamx

le problème c'est que si j'initialise les x0, y0 et z0 par des nombres qui différent de la solution par 0.01 par exemple alors la solution diverge

Salut
C'est probablement dû au problématique de l'imprécision des flottants dont j'ai parlé ici http://www.developpez.net/forums/d10...ontenant-float

Envoyé par hichamx

en plus chaque fois si je veux afficher le nombre d’itérations (compt) il m'affiche toujours 1
voici le code

Code c :

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
 
compt=1;
 
while( abs(x-x0)>epsilon  &&  abs(y-y0)>epsilon  &&  abs(z-z0)>epsilon )
 
{		
        x0=x;
        y0=y;
        z0=z;
 
        z=tab2[2][3] / tab2[2][2];
		y=((tab2[1][3]-(tab2[1][2]*z)))/tab2[1][1];
		x=((tab2[0][3]-(tab2[0][2]*z)-(tab2[0][1]*y)))/tab2[0][0];
 
        compt=compt+1;
 
}

Parce que la boucle ne fait qu'une seule itération. Probablement parce que tu n'as pas appliqué les lois de "De Morgan".
Je peux me tromper car je ne suis pas assez calé en maths pour vérifier ta façon de résoudre mais il me semble que tu dois quitter la boucle quand x correspond et aussi quand y correspond et aussi quand z correspond.
Donc, il s'ensuit que tu dois, inversement, rester dans la boucle tant que x ne correspond pas ou bien tant que y ne correspond pas ou bien quand z ne correspond pas. Déduction logique que De Morgan a si bien codifié avec son équation "non (a et b) = non(a) ou non(b)"

Code c :

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

while( abs(x-x0)>epsilon  ||  abs(y-y0)>epsilon  ||  abs(z-z0)>epsilon )

**hichamx** · 18/05/2011, 21h00

merci pour tes remarques mais le problème toujours existe: le programme ne fait qu'une seule itération

**diogene** · 19/05/2011, 17h04

Quelques remarques :
- De préférence, faire les calculs en double
-

Code :

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

while( abs(x-x0)>epsilon  ||  abs(y-y0)>epsilon  ||  abs(z-z0)>epsilon )

abs() s'applique à des entiers. Il convient ici d'utiliser fabs()

- Il est probable que ces formules sont fausses :

Code :

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

1
2
3
		z=tab2[2][3] / tab2[2][2];
		y=((tab2[1][3]-(tab2[1][2]*z)))/tab2[1][1];
		x=((tab2[0][3]-(tab2[0][2]*z)-(tab2[0][1]*y)))/tab2[0][0];

Quand on calcule y, on utilise la nouvelle valeur de z.
Quand on calcule x, on utilise les nouvelles valeurs de y et z.
C'est suspect.
Ceci étant, je n'ai rien compris à ces manipulations sur ces tableaux et le rapport avec l'algo de Newton - Raphson

- C'est sûr que l'algorithme est faux.
On calcule tab[] et tab2[] une bonne fois pour toutes avant de boucler. Les valeurs dans ces tableaux doivent dépendre du point courant, pas seulement du point initial.

**hichamx** · 24/05/2011, 18h33

D'ailleurs merci infiniment
j'ai bien modifié le programme, et si j'utilise les doubles et/ou la loi De Morgan il ne fait qu'une seule itération mais si je laisse les "floats" il fait n itérations ça veut dire il ne s’arrête pas, je suppose qu'il ne vérifie pas la condition d’arrêt

voici mon nouveau code

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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
 
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
float f (float x,float y,float z)
{ return (x*x+y*y+z*z-9);}
float g (float x,float y,float z)
{ return (x*x+y*y-3*x);}
float h (float x,float y,float z)
{ return (y+z-2);}
float df_1(float x,float y,float z)
{ return (2*x);}
float df_2(float x,float y,float z)
{ return (2*y);}
float df_3(float x,float y,float z)
{ return (2*z);}
float dg_1(float x,float y,float z)
{ return (2*x-3);}
float dg_2(float x,float y,float z)
{ return (2*y);}
float dg_3(float x,float y,float z)
{ return (0);}
float dh_1(float x,float y,float z)
{ return (0);}
float dh_2(float x,float y,float z)
{ return (1);}
float dh_3(float x,float y,float z)
{ return (1);}
 
 
 
 
 
int main()
{
 
    float x0,y0,z0,x[19],x1[19],x2[19],a[19][19], b[19],p,s,eps=0.001;
    int i,j,k,iter=1,n=3;
 
 
    printf("donner la valeur de x0 ");
scanf("%f",&x0);
 
printf("donner la valeur de y0 ");
scanf("%f",&y0);
 
printf("donner la valeur de z0 ");
scanf("%f",&z0);
 
do
{
 
    a[0][0]= df_1( x0,y0,z0);
    a[0][1]= df_2( x0,y0,z0);
    a[0][2]= df_3( x0,y0,z0);
    b[0]= (df_1( x0,y0,z0))*x0+(df_2( x0,y0,z0))*y0+(df_3( x0,y0,z0))*z0-f ( x0,y0,z0);
    a[1][0]= dg_1( x0,y0,z0);
    a[1][1]= dg_2( x0,y0,z0);
    a[1][2]= dg_3( x0,y0,z0);
    b[1]= (dg_1( x0,y0,z0))*x0+(dg_2( x0,y0,z0))*y0+(dg_3( x0,y0,z0))*z0-g ( x0,y0,z0);
    a[2][0]= dh_1( x0,y0,z0);
    a[2][1]= dh_2( x0,y0,z0);
    a[2][2]= dh_3( x0,y0,z0);
    b[2]= (dh_1( x0,y0,z0))*x0+(dh_2( x0,y0,z0))*y0+(dh_3( x0,y0,z0))*z0-h ( x0,y0,z0);
 
printf("%f",a[2][2]);
 
 
 
 
for(k=0;k<n-1;k++)
{
if (a[k][k]==0)
{
printf("\n\n * Un pivot nul ! => methode de Gauss non applicable\n\n");
system("PAUSE");main();
}
//réduction
for(i=k+1;i<n;i++)
{
p=a[i][k]/a[k][k];
for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
b[i]=b[i]-p*b[k];
}
}
//Résolution
for(i=n-1;i>=0;i--)
{
s=0;
for(j=i+1;j<n;j++)s=s+a[i][j]*x[j];
x[i]=(b[i]-s)/a[i][i];
}
 
for(i=0;i<n;i++)
{
for (j=0;j<n;j++) if (fabs(a[i][j])<eps) a[i][j]=0;
if (fabs(b[i])<eps) b[i]=0;
}
 
 
printf("\n-------------- Gauss -------------\n");
printf("\n * La matrice reduite :");
 
printf("\n\n");
for (i=0;i<n;i++)
{
printf(" [");
for (j=0;j<n;j++)
{
printf(" %.4f ",a[i][j]);
}
printf("] [ %.4f ]",b[i]);
printf("\n");
}
 
 
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 
iter++;
 
x0=x[0];
y0=x[1];
z0=x[2];
 
 
 
 
 
}while( (fabs(x[0]-x0)<eps)  &&  (fabs(x[1]-y0)<eps)  &&  (fabs(x[2]-z0)<eps) ) ;
printf("le nombre d'iterations est %d\n",iter);
 
system("PAUSE");
}

**Sve@r** · 03/06/2011, 14h06

Envoyé par hichamx

D'ailleurs merci infiniment
j'ai bien modifié le programme, et si j'utilise les doubles et/ou la loi De Morgan il ne fait qu'une seule itération mais si je laisse les "floats" il fait n itérations ça veut dire il ne s’arrête pas, je suppose qu'il ne vérifie pas la condition d’arrêt

Salut

A ce niveau là, seul un printf() de tes valeurs dans la boucle assortie éventuellement d'un getchar() pour te laisser le temps de voir le résultat permettra de lever le doute.

Mais les équations de "De Morgan" ne sont pas une formule magique et quelque peu aléatoire mais des lois mathématiques clairement établies et 100% efficaces. Le contraire de "avoir A ET B" sera "ne pas avoir A OU ne pas avoir B". Et le contraire de "avoir A OU B" sera "ne pas avoir A ET ne pas avoir B".
Exemple:
je veux une personne majeure et en âge de travailler => while (age >= 18 et age < 60)
je veux l'inverse => while (age < 18 ou age >= 60)
Si tu veux t'en convaincre, une table de vérité avec A, B, NON(A ET B), NON(A) OU NON(B) te montrera l'égalité des 2 dernières colonnes. Et idem avec A, B, NON(A OU B), NON(A) ET NON(B).

Donc toi tu as écrit while( (fabs(x[0]-x0)<eps) && (fabs(x[1]-y0)<eps) && (fabs(x[2]-z0)<eps) )
Ce qui signifie "je reste dans la boucle tant que x[0] quasi égal à x0 et x[1] quasi égal à y0 et x[2] quasi égal à z0"
Donc, inversement, tu quittes la boucle si "x[0] différent de x0 ou x[1] différent de y0 ou x[2] différent de z0". Comme je te l'ai dit, mon niveau mathématique n'est pas assez fort pour estimer si cette condition est mathématiquement correcte mais mon intuition me dit que non...

PS: tu devrais séparer tes tâches en fonctions. Une fonction qui saisit, une qui calcule. Ca améliorera ton analyse et ton débug...

**souviron34** · 03/06/2011, 14h38

Envoyé par hichamx

D'ailleurs merci infiniment
j'ai bien modifié le programme, et si j'utilise les doubles et/ou la loi De Morgan il ne fait qu'une seule itération mais si je laisse les "floats" il fait n itérations ça veut dire il ne s’arrête pas, je suppose qu'il ne vérifie pas la condition d’arrêt

le passage de paramètres en C est en double lorsque le nombre est à virgule flottante, quel que soit le type d'entrée..

Si le type d'entrée est un float, c'est converti en double, puis reconverti. Même chose en sortie.. Il y a donc perte.

Pour tout calcul passant par des fonctions, il est plus sûr d'utiliser des doubles, d'autant plus que c'est le type naturel des opérations de maths (sqrt, pow, log, exp).

**diogene** · 03/06/2011, 19h10

le passage de paramètres en C est en double lorsque le nombre est à virgule flottante, quel que soit le type d'entrée..

Non, rien n'oblige à cela, sauf dans le cas où le type du paramètre n'est pas spécifié par le prototype.
Les opérations sur float peuvent être faites en double (sauf cast et assignation), mais ce n'est pas une obligation.

Si le type d'entrée est un float, c'est converti en double, puis reconverti. Même chose en sortie.. Il y a donc perte.

Non : Dans la conversion double->float, si la valeur peut être représentée par un float, la conversion ne change pas la valeur et c'est le cas lors de la suite de conversion float->double->float puisque la conversion float->double n'altère pas la valeur.