Intégrer une équation

**lisenette** · 07/04/2014, 11h51

Bonjour,
je voudrais intégrer une équation deux fois sur Matlab ; je sais qu'il existe une fonction "ODE 45" mais je ne sais pas comment l'utiliser.

Invité · 07/04/2014, 12h03

Bonjour,

Tu trouveras sa documentation en tapant doc ode45 dans MATLAB.
N'hésite pas non plus à effectuer une recherche sur le forum : l'utilisation de cette fonction a été abordée de nombreuses fois.

**lisenette** · 07/04/2014, 14h17

j'essaye d'utiliser la syntaxe donnée sur matlab

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
 
xf0=0.2;
yf0=0.2;
vxf0=0.791;
vyf0=0.791;
xi=xf0;
yi=yf0;
vxi=vxf0;
vyi=vyf0;
T=0;
V0=[xi,yi,vxi,vyi];
 
ll=T+Tsup;
    %pendant la phase de double-support
   [TT,V]=ode45(@sysmodel,[T:0.02:ll],V0);

et voici l'erreur que sa me donne

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
??? Error using ==> sysmodel
Too many input arguments.

Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> ode at 51
   [TT,V]=ode45(@sysmodel,[T:0.02:ll],V0);

la fonction qui contient le modèle à intégrer est:

Code :

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

1
2
3
4
5
6
7
8
9
 
%modele dynamique
function [xpp,ypp]=sysmodel()
%global pxmp  pymp
g=9.8;
z=0.8;
xpp=(g/z)*(x);
ypp=(g/z)*(y);
end

Invité · 07/04/2014, 14h43

La fonction dont le handle est passé à ode45 : sysmodel dans ton cas doit posséder deux arguments en entrée, et un en sortie. Soit de la forme (pour reprendre les noms de la doc) :

Code :

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

function dy = sysmodel(t,y)

dy est le vecteur colonne de sortie, exprimé en fonction des entrées t et y selon le système

y₁' = f(t, y₁, y₂, ...)
y₂' = f(t, y₁, y₂, ...)
...

**lisenette** · 07/04/2014, 15h53

c'est toujours la même erreur si j'utilise votre syntaxe sachant que la fonction n'admet pas d'argument et que x et y sont des vecteurs quelconques.

Code :

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

1
2
3
4
5
6
7
8
9
10
 
%modele dynamique
function acc=sysmodel()
%global pxmp  pymp
g=9.8;
z=0.8;
xpp=(g/z)*(x);
ypp=(g/z)*(y);
acc=[xpp,ypp];
end

Invité · 07/04/2014, 16h11

Envoyé par Winjerome

doit posséder deux arguments en entrée, et un en sortie.

Code :

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

function dy = sysmodel(t,y)

Des entrées qu'il faudra utiliser d'ailleurs

À la place de x et y j'imagine...

N'hésite pas à exécuter les exemples de la documentation au pas à pas pour bien comprendre le fonctionnement.

**lisenette** · 08/04/2014, 09h35

merci pour vos réponses.

je voudrais savoir si ODE45 peut me donner à la fois la solution en vitesse et en position ?

**lisenette** · 08/04/2014, 10h25

voici donc mon programme complet:

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
 
% clc,
% clear all,
% close all, 
%declaration des constantes
g=9.8;
Tsup=0.8;%temps de support
z=0.8;%hauteur du centre de masse considérée constante pour un modele LIP
a=10;
b=1;
Tc=sqrt(z/g);
Tdbl=0.1;
pas=0.02;
sx0=0;
sy0=0.2;
sx=[0.0 0.3 0.3 0.3 0];%longueur des pas
sy=[0.2 0.2 0.2 0.2 0.2];%largeur des pas
 
%CI1=[-0.2,0.791];
pxmp=0;
pymp=0;
px=0;
py=0;
pxprevious=0;
pyprevious=0;
n=length(sx);
L=0.06;
%etape2:
%c'est valeures sont choisi apres calul ou j'ai remplacé t=0  et r=0
xf0=0.2;
yf0=0.2;
vxf0=0.791;
vyf0=0.791;
xi=xf0;
yi=yf0;
vxi=vxf0;
vyi=vyf0;
i=0;
T=0;
 
for r=1:length(sx)
    t=T+Tsup;
    V0x=[xi,vxi];
    V0y=[yi,vyi];
   [TT,Vx]=ode45(@sysmodelx,[T:0.02:t],V0x);
   [SS,Vy]=ode45(@sysmodely,[T:0.02:t],V0y);
 
%Etape5:
        px=pxprevious+sx0;
        py=pyprevious-((-1)^r)*sy0;%si on choisit le pied gauche comme point de départ de la marche on remplace -(-1)^n par +(-1)^n
        pxprevious=px;
        pyprevious=py;
 
%         hold on,
% plot(px,py,'r+'),grid on; 
%Etape6:les caracteristiques de la primitive de marche suivante 
%Ajouter une condition
 
        xbar=sx0/2;
        ybar=((-1)^r)*sy0/2;
        vxbar=(cosh(Tsup/Tc)+1)/(Tc*sinh(Tsup/Tc))*xbar;
        vybar=(cosh(Tsup/Tc)-1)/(Tc*sinh(Tsup/Tc))*ybar;
 
%Etape7:calcul de l'etat cible
        xd=px+xbar;
        vxd=vxbar;
        yd=py+ybar;
        vyd=vybar;
        C=cosh(Tsup/Tc);
        S=sinh(Tsup/Tc);
        D=a*(C)^2+b*(S/Tc)^2;     
 
%Etape8:calcul de la position modifiée du pied
pxm=(a*C*(xd-S*Tc*vxd)+b*S/Tc*(vxd-C*vxi))/D;
pym=(a*C*(yd-S*Tc*vyd)+b*S/Tc*(vyd-C*vyi))/D;
pxmp=pxm;
pymp=pym;
T=ll;
sx0=sx(r+1);
sy0=sy(r+1);    
xi=X;
vxi=Vx;
yi=Y;
vyi=Vy; 
 
 
end

ainsi que es fonctions à utiliser dans ODE45 :

Code :

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

1
2
3
4
5
6
7
 
%modele dynamique suivant l'axe x
function xpp=sysmodelx(t,x)
global g  z pxm
xpp=(g/z)*(x-pxm);
 
end

Code :

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

1
2
3
4
5
6
7
 
%modele dynamique suivant l'axe x
function ypp=sysmodely(t,y)
global g  z pym
ypp=(g/z)*(y-pym);
 
end

mais j'ai toujours une erreur :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
??? Error using ==> minus
Matrix dimensions must agree.

Error in ==> sysmodelx at 4
xpp=(g/z)*(x-pxm);

Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> ode at 44
   [TT,Vx]=ode45(@sysmodelx,[T:0.02:t],V0x);

la variable de temps n’apparaît pas dans mon équation,comment substituer cette variable ?
mon système est de la forme : xpp=a*(x-xc) et ypp=a*(y-yc)
avec xpp et ypp les dérivées secondes de x et y respectivement.

éclairez-moi je vous pris.

Invité · 08/04/2014, 12h54

Ton équation étant d'ordre 2, il faut faire intervenir la dérivée d'ordre 1 pour se retrouver avec deux équations d'ordre 1. Le temps n'intervenant pas dans ton équation, l'entrée t ne retrouvera non-utilisée ici.
Tu peux t'inspirer de cette discussion.
[EDIT]Voir ce message

**lisenette** · 09/04/2014, 11h24

sa ne marche toujours pas et franchement je n'y comprends rien.
bon, j'ai mis le système sous forme d’Équation d’état mais c'est la même erreur:
mon programme principal est:

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
 
% clc,
% clear all,
% close all, 
%declaration des constantes
g=9.8;
Tsup=0.8;%temps de support
z=0.8;%hauteur du centre de masse considérée constante pour un modele LIP
a=10;
b=1;
Tc=sqrt(z/g);
Tdbl=0.1;
pas=0.02;
sx0=0;
sy0=0.2;
sx=[0.0 0.3 0.3 0.3 0];%longueur des pas
sy=[0.2 0.2 0.2 0.2 0.2];%largeur des pas
 
%CI1=[-0.2,0.791];
pxmp=0;
pymp=0;
px=0;
py=0;
pxprevious=0;
pyprevious=0;
n=length(sx);
L=0.06;
%etape2:
%c'est valeures sont choisi apres calul ou j'ai remplacé t=0  et r=0
xf0=0.2;
yf0=0.2;
vxf0=0.791;
vyf0=0.791;
xi=xf0;
yi=yf0;
vxi=vxf0;
vyi=vyf0;
i=0;
T=0;
    X=[xi,vxi];
    Y=[yi,vyi];
    X0=X;
    Y0=Y;
 
for r=1:length(sx)
    t=T+Tsup;
    X=[xi,vxi];
    Y=[yi,vyi];
   [TT,x]=ode45(@sysmodelx,[T t],X0);
   [SS,y]=ode45(@sysmodely,[T t],Y0);
 
%Etape5:
        px=pxprevious+sx0;
        py=pyprevious-((-1)^r)*sy0;%si on choisit le pied gauche comme point de départ de la marche on remplace -(-1)^n par +(-1)^n
        pxprevious=px;
        pyprevious=py;
 
%         hold on,
% plot(px,py,'r+'),grid on; 
%Etape6:les caracteristiques de la primitive de marche suivante 
%Ajouter une condition
 
        xbar=sx0/2;
        plot(px,py);
        plot(pmx,pmy);
        ybar=((-1)^r)*sy0/2;
        vxbar=(cosh(Tsup/Tc)+1)/(Tc*sinh(Tsup/Tc))*xbar;
        vybar=(cosh(Tsup/Tc)-1)/(Tc*sinh(Tsup/Tc))*ybar;
 
%Etape7:calcul de l'etat cible
        xd=px+xbar;
        vxd=vxbar;
        yd=py+ybar;
        vyd=vybar;
        C=cosh(Tsup/Tc);
        S=sinh(Tsup/Tc);
        D=a*(C)^2+b*(S/Tc)^2;     
 
%Etape8:calcul de la position modifiée du pied
pxm=(a*C*(xd-S*Tc*vxd)+b*S/Tc*(vxd-C*vxi))/D;
pym=(a*C*(yd-S*Tc*vyd)+b*S/Tc*(vyd-C*vyi))/D;
pxmp=pxm;
pymp=pym;
T=T+Tsup;
sx0=sx(r+1);
sy0=sy(r+1);    
xi=X(1);
vxi=X(2);
yi=Y(1);
vyi=Y(2); 
 
 
end

et pour les fonctions:

Code :

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

1
2
3
4
5
6
7
8
9
 
%modele dynamique suivant l'axe x
function xpp=sysmodelx(x,X)
global g  z pxm
xpp=zeros(2,1);
xpp(1)=X(2);
xpp(2)=(g/z)*(X(1)-pxm);
 
end

Code :

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

1
2
3
4
5
6
7
8
9
 
%modele dynamique suivant l'axe x
function ypp=sysmodely(y,Y)
global g  z pym
ypp=zeros(2,1);
ypp(1)=Y(2);
ypp(2)=(g/z)*(Y(1)-pym);
 
end

Invité · 09/04/2014, 11h42

Pour utiliser des variables globales, il faut les déclarer globales partout où elles sont utilisées, y compris dans ton programme principal.
même cela fait, tu ne fixes la valeur de pxm qu'une fois les lignes ode passées, cette variable reste donc vide, ce qui pose évidemment problème.

N'hésite pas à déboguer ton programme, mettre des breakpoints là où il y a des problèmes pour regarder leur différentes valeurs afin de comprendre le pourquoi du comment.

**lisenette** · 09/04/2014, 17h52

maintenant je n'ai plus d'erreur au niveau de ode45 merci .
par contre en dessinant mes trajectoires à savoir x en fonction de y je devais obtenir des sinusoïdes mais moi j'obtiens des droites je me demande si se n'est pas une erreur de logique sur ode45 ?

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
 
clc,
clear all,
close all, 
%declaration des constantes
global g z pxm pym
g=9.8;
Tsup=0.8;%temps de support
z=0.8;%hauteur du centre de masse considérée constante pour un modele LIP
a=10;
b=1;
Tc=sqrt(z/g);
Tdbl=0.1;
pas=0.02;
sx0=0;
sy0=0.2;
sx=[0.0 0.3 0.3 0.3 0];%longueur des pas
sy=[0.2 0.2 0.2 0.2 0.2];%largeur des pas
sx(6)=0;
sy(6)=0.2;
 
%CI1=[-0.2,0.791];
pxmp=0;
pymp=0;
pxm=0;
pym=0;
pxprevious=0;
pyprevious=0;
n=length(sx);
L=0.06;
%etape2:
%c'est valeures sont choisi apres calul ou j'ai remplacé t=0  et r=0
xf0=0.2;
yf0=0.2;
vxf0=0.791;
vyf0=0.791;
xi=xf0;
yi=yf0;
vxi=vxf0;
vyi=vyf0;
i=0;
T=0;
    X=[xi,vxi];
    Y=[yi,vyi];
    X0=X;
    Y0=Y;
 
for r=1:5
    for t=T+Tsup;
    X=[xi,vxi];
    Y=[yi,vyi];
   [TT,x]=ode45(@sysmodelx,[T t],X0);
   [SS,y]=ode45(@sysmodely,[T t],Y0);
   hold on,
   subplot(1,2,1),plot(x(:,1),y(:,1));
   %plot(SS,y(:,1))
   %,title('x-distance')
   %subplot(1,2,2),plot(TT,x(:,2)),plot(SS,y(:,2),title('x-velocity');
 
%Etape5:
        px=pxprevious+sx0;
        py=pyprevious-((-1)^r)*sy0;%si on choisit le pied gauche comme point de départ de la marche on remplace -(-1)^n par +(-1)^n
        pxprevious=px;
        pyprevious=py;
 
%         hold on,
% plot(px,py,'r+'),grid on; 
%Etape6:les caracteristiques de la primitive de marche suivante 
%Ajouter une condition
 
        xbar=sx0/2;
        plot(px,py);
        %plot(pmx,pmy);
        ybar=((-1)^r)*sy0/2;
        vxbar=(cosh(Tsup/Tc)+1)/(Tc*sinh(Tsup/Tc))*xbar;
        vybar=(cosh(Tsup/Tc)-1)/(Tc*sinh(Tsup/Tc))*ybar;
 
%Etape7:calcul de l'etat cible
        xd=px+xbar;
        vxd=vxbar;
        yd=py+ybar;
        vyd=vybar;
        C=cosh(Tsup/Tc);
        S=sinh(Tsup/Tc);
        D=a*(C)^2+b*(S/Tc)^2;     
 
%Etape8:calcul de la position modifiée du pied
pxm=(a*C*(xd-S*Tc*vxd)+b*S/Tc*(vxd-C*vxi))/D;
pym=(a*C*(yd-S*Tc*vyd)+b*S/Tc*(vyd-C*vyi))/D;
T=T+Tsup;
sx0=sx(r+1);
sy0=sy(r+1);    
xi=X(1);
vxi=X(2);
yi=Y(1);
vyi=Y(2); 
 
    end    
end 
hold off,

Code :

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

1
2
3
4
5
6
7
8
9
 
%modele dynamique suivant l'axe x
function xpp=sysmodelx(x,X)
global g  z pxm
xpp=zeros(2,1);
xpp(1)=X(2);
xpp(2)=(g/z)*(X(1)-pxm);
 
end

Code :

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

1
2
3
4
5
6
7
8
9
 
%modele dynamique suivant l'axe x
function ypp=sysmodely(y,Y)
global g  z pym
ypp=zeros(2,1);
ypp(1)=Y(2);
ypp(2)=(g/z)*(Y(1)-pym);
 
end

Invité · 09/04/2014, 18h12

J'ai un doute sur les lignes qui précèdent :

Code :

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

1
2
        X=[xi,vxi];
        Y=[yi,vyi];

n'est serait-ce pas X0 et Y0 ?

Pense aussi à rajouter les abscisses :

Code :

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

plot(TT, x(:,1), SS, y(:,1));

Sinon la logique me semble bonne.

**lisenette** · 10/04/2014, 11h48

je ne sais pas pourquoi j'ai toujours des trajectoires bizarres et non polynomiale enfin curviligne ?
je bugg

voici le trajectoire que j'obtiens !!!!! Nom : trajectoire.png
Affichages : 229
Taille : 44,4 Ko

Invité · 10/04/2014, 12h36

C'est à dire continues ? Si oui, ce sont les conditions initiales X0, Y0 que tu dois revoir.

**lisenette** · 10/04/2014, 14h26

merci winjerome.
je voudrais poser une question concernant les conditions initiales:
pourquoi influent-elle sur la trajectoire et comment les choisir (critère).

Invité · 10/04/2014, 14h42

Les conditions initiales correspondent à la position et la vitesse que tu donnes à ton objet au départ.
Imagine une balle que tu lances vers le haut ou vers le bas, elle n'aura pas du tout la même trajectoire dans un cas comme dans l'autre

Dans ton cas, tu sembles effectuer un découpage selon plusieurs intervalles de temps qui se suivent. Si tu redonnes exactement les mêmes positions et vitesses à chaque fois, ne sois pas étonné que ton objet reprenne exactement la même trajectoire.
Au contraire si tu donnes les dernières valeurs prises pour le "prochain lancé", là ton objet continuera sur sa lancée.

**lisenette** · 10/04/2014, 16h21

je viens de déterminer mes conditions initiales me permettant d'avoir l'allure voulue mais je l'ai fait par tâtonnement, y a-t-il une méthode pour les déterminer avant de dessiner.

merci encore pour ton aide winjerome.

Invité · 10/04/2014, 17h02

Envoyé par Winjerome

Au contraire si tu donnes les dernières valeurs prises pour le "prochain lancé", là ton objet continuera sur sa lancée.

Code :

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

1
2
3
4
5
6
7
8
9
10
X0 = ...
Y0 = ...
 
for ...
    [TT,x]=ode45(@sysmodelx,[T t],X0);
    [SS,y]=ode45(@sysmodely,[T t],Y0);
 
    X0 = x(end,:);
    Y0 = y(end,:);
end

Intégrer une équation

MATLAB

Discussions similaires

Partager

Partager