Jacobian dans une fonction

Bonjour,

Je dois faire une fonction qui me calcul le jacobien, et qui me donne et sortie la matrice sous cette forme bien précise :

Code:

---

J=[0,1;-2000*x(1)*x(2)-1,1000-1000*x(1)^2];

---

Je calcul le jacbien de la manière suivante :

Code:

---

1 2 3 4

syms x1 x2 h=[x2,1000*(1-x1^2)*x2-x1]; v=[x1 x2]; J=jacobian(h,v);

---

Le problème c'est que j'obtiens le jacobien avec des x1 et x2 au lieu de l'avoir sous forme vectorielle avec x(1) et x(2). Ca paraît bête mais j'ai essayé beaucoup de truc, notamement avec "subs" mais rien n'y fait...

Est-ce que quelqu'un saurait comment à la fois calculer le jacobien dans ma fonction, et avoir la sortie sous la bonne forme?

salut

tu peux simplifier ton écriture comme ceci

Code:

---

v = sym('x',[1,2]);

---

mais ca change rien, il te donnera tjrs le résultat en fonction de x1 et x2 et jamais en fonction de v(1) et v(2)

par contre tu peux ajouter

Code:

---

1 2	J = subs(J,'x1','x(1)'); J = subs(J,'x2','x(2)');

---

Fabien

Alors j'ai essayé :

Code:

---

1 2	J = subs(J,'x1','x(1)'); J = subs(J,'x2','x(2)');

---

Mais j'obitens ça comme message :

"Conversion to logical from sym is not possible."

:(

étrange :

Code:

---

1 2 3 4 5 6 7 8 9

>> syms x1 x2 h=[x2,1000*(1-x1^2)*x2-x1]; v=[x1 x2]; J=jacobian(h,v)   J =   [ 0, 1] [ - 2000*x1*x2 - 1, 1000 - 1000*x1^2]

---

Code:

---

1 2 3 4 5 6 7 8

>> J = subs(J,'x1','x(1)'); J = subs(J,'x2','x(2)')   J =   [ 0, 1] [ - 2000*x(1)*x(2) - 1, 1000 - 1000*x(1)^2]

---

Ahh ok, alors maintenant que j'essai dans le terminal matlab, ça marche aussi.

Mais le problème c'est que dans le reste de mon script ça ne marche pas...

Code:

---

1 2 3 4 5 6 7 8 9 10 11

Voilà la fonction complète :   function J=FJac(t, x)   syms x1 x2 h=[x2,1000*(1-x1^2)*x2-x1]; v=[x1 x2]; J=jacobian(h,v);   J = subs(J,'x1','x(1)'); J = subs(J,'x2','x(2)');

---

et quand je tape par exemple dans le terminal :

Code:

---

1 2	vect=[1 1] FJac(1,vect)

---

j'ai :

Code:

---

1 2 3 4

ans =   [ 0, 1] [ - 2000*x(1)*x(2) - 1, 1000 - 1000*x(1)^2]

---

j'en déduit qu'il ne reconnait pas mon vecteur "vect" comme argument de ma fonction, ou quelque chose comme ça puisque la réponse devrait être :

Code:

---

1 2 3 4

ans =   0 1 -2001 0

---

Je trouve ça assez bizarre...

salut

alors la "solution" que je t'ai donnée n'est pas terrible car elle fait grosso modo une substitution de string :
je pensais que tu voulais faire de l'affichage !!

pour faire le calcul, il faut l'écrire comme ca :
1/ le calcul en symbolique

Code:

---

1 2 3 4 5 6 7 8 9

>> syms x1 x2 h=[x2,1000*(1-x1^2)*x2-x1]; v=[x1 x2]; J=jacobian(h,v)   J =   [ 0, 1] [ - 2000*x1*x2 - 1, 1000 - 1000*x1^2]

---

2/le calcul avec les valeurs réelles :

Code:

---

1 2 3 4 5 6 7

>> vect = [1 1]; >> subs(J,v,vect)   ans =   0 1 -2001 0

---

Fabien

Merci beacoup de ton aide, mais malheureusement ça ne fonctionne pas encore... en fait, ce que j'aurai du dire dès début, c'est que je calcule le jacobian dans le contexte de la résolution d'une équation différentielle. Je n'ai donc pas directement accès au input de ma fonction qui calcule le Jacobian. Voilà les 3 scripts qui constitue mon "programme"

1) l'énoncé de mon système d'équations différentielles

Code:

---

1 2 3 4 5 6

function f = vdp_custom(t, x)   f(1) = x(2); f(2) = 1000*(1 - x(1)^2)*x(2) - x(1);   f = f';

---

2)le calcul du jacobian (qui ne fonctionne pas dans ce contexte)

Code:

---

1 2 3 4 5 6 7 8 9 10

function J=FJac(t, x)   syms x1 x2 h=[x2,1000*(1-x1^2)*x2-x1]; v=[x1 x2]; J=jacobian(h,v);     J = subs(J,'x1','x(1)'); J = subs(J,'x2','x(2)');

---

3)la résolution de l'équa diff

Code:

---

1 2 3 4 5 6 7

t_init=0; t_final=3000; x1_0=2; x2_0=1;   options_struct = odeset('RelTol', 1E-9, 'AbsTol', 1E-6,'Stats','on','Jacobian',@FJac); [t, x] = ode23s(@vdp_custom,[t_init t_final],[x1_0 x2_0], options_struct);

---

On voit ici que les inputs de Fjac sont les valeurs de t et x qui sont calculés par incrémentations successives dans le solver ode23s, auquel je n'ai pas accès en cours de résolution...

Il faudrait donc que je puisse faire une substitution directement dans ma fontion Fjac, et que le programme reconnaisse mes x(1) et x(2) comme arguements de Fjac.

Au pire des cas, j'écris direment mon jacobian dans Fjac et c'est réglé... mais j'aurai aimé le calculer dans la fonction pour pouvoir utilisé ma fonction pour d'autres problèmes.

et si tu écris ta fonction FJac comme ceci :

Code:

---

1 2 3 4 5 6 7 8 9

function J=FJac(t, x)   syms x1 x2 h=[x2,1000*(1-x1^2)*x2-x1]; v=[x1 x2]; J=jacobian(h,v);     J = subs(J,v,x);

---

Non... maintenant j'ai ça comme erreur :

Code:

---

OLD and NEW must have the same number of entries.

---

Je pense que comme il ne sait pas encore quelle taille feront les arguments de Fjac il bloque... quel merde

salut

je sais pas comment tu t'y prend, mais ça marche chez moi :

le code :

Code:

---

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

function test   t_init=0; t_final=3000; x1_0=2; x2_0=1;   options_struct = odeset('RelTol', 1E-9, 'AbsTol', 1E-6,'Stats','on','Jacobian',@FJac); [t, x] = ode23s(@vdp_custom,[t_init t_final],[x1_0 x2_0], options_struct);   whos t x   figure plotyy(t,x(:,1),t,x(:,2))   function f = vdp_custom(t, x)   f(1) = x(2); f(2) = 1000*(1 - x(1)^2)*x(2) - x(1);   f = f';   function J=FJac(t, x)   syms x1 x2 h=[x2,1000*(1-x1^2)*x2-x1]; v=[x1 x2]; J=jacobian(h,v);     J = subs(J,v,x);

---

l'appel et le résultat

Code:

---

1 2 3 4 5 6 7 8 9 10 11

>> test 10577 successful steps 8 failed attempts 31749 function evaluations 10577 partial derivatives 10585 LU decompositions 31755 solutions of linear systems Name Size Bytes Class Attributes   t 10578x1 84624 double x 10578x2 169248 double

---

fabien

Ah oui?

Alors là j'avoue que je comprends pas... J'ai copié-collé ton code et essayé sur 2 ordinateurs différents et j'ai sur les deux le même message d'erreur...

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Error using sym/subs>normalize (line 201) OLD and NEW must have the same number of entries.   Error in sym/subs>mupadsubs (line 136) [X2,Y2,symX,symY] = normalize(X,Y); %#ok   Error in sym/subs (line 124) G = mupadsubs(F,X,Y);   Error in test2>FJac (line 31) J = subs(J,v,x); Error in ode23s (line 221) dfdy = feval(Jac,t,y,Jargs{:});   Error in test2 (line 9) [t, x] = ode23s(@vdp_custom,[t_init t_final],[x1_0 x2_0], options_struct);

---

Peut être que c'est lié à la version de matlab... Merci du temps que tu as consacré pour m'aider en tout cas!

salut

je suis sous 2011b

sinon, met un point d'arrêt sur la ligne en cause (le subs) et regarde les tailles et valeurs (nan?) de x
peut être une évolution de ode23s entre ta version et la mienne, que peut être tu peux gérer en interceptant le problème

Fabien

Moi j'ai la 2013a. J'ai vu mon assistant aujourd'hui et il avait l'air de dire que c'était pas non plus possible sur notre version.

En tout cas merci beaucoup pour ton aide, je vais cloturer l'affaire maintenant ^^