Discussion :

[cpalperou]

Bonjour,
je dois résoudre une 'équation différentielle' du premier ordre très simple du type dy/dx = f(x) avec la fonction ode23 ou ode45 (c'est imposé)
J'ai une grille pour les valeurs de x sous la forme d'un vecteur: (x1,x2,...,xn) et les valeurs de f correspondantes: (f(x1),f(x2),...).
Quelqu'un peut il m'aider SVP?
merci

[duf42]

Bonjour,

Peux-tu nous montrer ce que tu as commencé à faire et les difficultés précises que tu rencontres?

Duf

[cpalperou]

Bonjour,
mon problème vient, je pense de la simplicité du problème.
Je comprends bien comment utiliser ode45 avec un système quelconque d'équa diff du genre dy1/dx=f(y1,y2,x) et dy2/dt=g(y1,y2x)
Mais là, j'ai juste dy/dx=f(x) donc, je fais un blocus total!

peut être ceci:

Code :

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

1
2
3
4
 
theta = [0 1 2 3 4 5]; % c'est juste un exemple, la réalité est plus cplx
sqq=@theta
[T,Y]=ode45(sqq,[0 100],12);

[cpalperou]

Bonjour,
je n'y arrive toujours pas et je n'ai pas de réponses de votre part. Je m'explique peut être mal?
Voici plus explicitment mon problème:

j'ai un vecteur Tr de dimension 500, un vecteur L de dimension 500. Chaque valeur de Tr et de L étant fonction du temps t. Je calcule la fonction

Code :

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

y=1./(L.*Tr) - 1./L;

Jusque là, pas de problème.
Je sais ensuite que cette expression est égale à une dérivée du temps, soit da/dt = y (avec da /dt la dérivée de la fonction a(t).

Je dois alors retrouver la fonction a(t) en utilisant un solveur d'équation différentielle (ode45 par exemple), ceci est imposé.

J'écris donc une fonction rigid.m ainsi:

Code :

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

1
2
3
 
function dy = rigid(t,y,Tr,L)
dy = 1./(L.*Tr)-1./L;

puis:

Code :

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

1
2
3
4
 
tspan=[0 20];
IC = 1e-2;
[t,y]=ode23(@(t,y) rigid(t,y,Tr,L),tspan,IC);

et j'obtiens une erreur:

Code :

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

1
2
3
4
5
Error using odearguments (line 93)
@(T,Y)RIGID(T,Y,Tr,L) returns a vector of length 500, but the length of initial conditions vector is 1. The vector returned by @(T,Y)RIGID(T,Y,Tr,L) and the initial conditions vector must have the same number of elements.
 
Error in ode23 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

SVP, merci

[Invité]

Bonjour,

Le rôle de la fonction odefun - ou rigid dans ton cas - est de retourner la valeur dy correspondant aux entrées (t,y) que lui donne la fonction odeXX en interne, pas tout le temps le même vecteur quelques soient ces entrées.

[cpalperou]

Merci winjérome,
je pense avoir bien compris le rôle de la fonction odefun (rigid dans mon cas). Je n'arrive toutefois pas à l'utiliser correctement!
j'ai changé ma fonction rigid.m ainsi:

Code :

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

1
2
3
function dy = rigid(t,y,Tr,L)
y = 1./(L.*Tr)-1./L;
dy = y(1);

je n'ai plus d'erreur mais le résultat de

Code :

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

1
2
3
4
 
tspan=[0 20];
IC = 1e-2;
[t,y]=ode23(@(t,y) rigid(t,y,Tr,L),tspan,IC);

me donne un vecteur y de 11 valeurs toutes égales à IC=1e-2.
J'm'en sors pas punaise!! J'y retourne, je dois être bête, c'est pas possible!!

[Invité]

Comme tu l'as dit

Chaque valeur de Tr et de L étant fonction du temps t.

Ce temps que justement tu as en entrée, donc sers t'en

[cpalperou]

Encore merci,
j'ai essayé ceci:

Code :

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

1
2
3
4
5
6
7
 
function dy = rigid(t,y,Tr,L)
y = 1./(L.*Tr(:,1))-1./L;
dy=zeros(numel(y),1);
for i=1:numel(y)
    dy(i,1) = y(i);
end

ainsi, j'ai, pour chaque valeur de t le dy correspondant mais quand je fais ode23,
j'ai le même problème de condition initial:

Code :

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

1
2
3
4
5
Error using odearguments (line 93)
@(T,Y)RIGID(T,Y,Tr,L) returns a vector of length 500, but the length of initial conditions vector is 1. The vector returned by @(T,Y)RIGID(T,Y,Tr,L) and the initial conditions vector must have the same number of elements.
 
Error in ode23 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Je continue de chercher

[Invité]

De 1), tu écrases le y d'entrée :

Code :

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

1
2
function dy = rigid(t,y,Tr,L)
y = 1./(L.*Tr(:,1))-1./L;

Et bien que tu ne te serves pas de cette entrée, autant garder les choses au clair
Ensuite, peux-tu me dire en quoi le code que tu me montres fait dépendre la sortie dy de l'unique valeur que contient l'entrée t ? Et en quoi cela diffère de ce que tu faisais dans ton 3ème message ?

[cpalperou]

J'avoue être perdu! Je n'arrive même pas à répondre à ta question! Je retourne y réfléchir

[cpalperou]

Bonjour,
je n'y arrive toujours pas, je suis complètement bloqué. Pouvez vous me débloquer svp?
J'en suis là:

Code :

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

1
2
function dy=rigid(t,y,Tr,L)
dy = 1./(L.*Tr)-1./L;

puis :

Code :

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

1
2
3
tspan=[0 20];
IC = 1e-2;
[t,y]=ode23(@(t,y) rigid(t,y,Tr,L),tspan,IC);

et: erreur !! RRRRrrrrRrrrr Je n'arrive pas à faire intervenir le temps t dans tout ça (comme me le conseille Winjérome!)

[Jerome Briot]

Il y a bien longtemps que je n'ai pas utilisé ces fonctions, mais ne faudrait-il pas s'appuyer sur l'exemple n°3 de la documentation ?

[Invité]

Tout à fait, suffisait de lire la documentation

[tachmou]

Salut !

Oulah, c'est un joyeux bordel ici... bon, reprenons les choses calmement.

Tu dois utiliser ode45 disons. Très bien, ouvre la doc des solver à côté pour l'avoir en parallèle à mes explications. Il y a des choses basiques à comprendre et à maîtriser.

Code :

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

[T,Y] = solver(odefun,tspan,y0)

Le solver, ici, sera ode45. Donc tu vas utiliser quelque chose du style :

Code :

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

[T,Y] = ode45(odefun,tspan,y0)

Bien. Maintenant, regardons le tableau en dessous (dans la doc)... Tu as 3 paramètres à spécifier en entrée de ode45 :

- y0 : conditions initiales de ton vecteur solution. Si tu as un seul paramètre, c'est donc une simple constante
- tspan : quelque chose du style [t0,tf] par exemple
- odefun : là c'est plus compliqué... et c'est là que tu butes.

odefun doit être une fonction qui doit être de la forme :

Code :

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

dydt = odefun(t,y)

>>> Il est impossible de mettre autre chose que t et y en entrée de cette fonction ! <<< Et la seule sortie possible est le résultat du calcul de la dérivée de y par rapport au temps (dydt) !

Donc, dydt est l'expression de y'(t) en fonction de y(t) et de t !

Conclusion...

Première étape : donne nous l'expression de y'(t) = f( t , y(t) ) - qui doit être valable quel que soit t.

Deuxième étape : écris une fonction de la forme :

Code :

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

1
2
3
4
5
dydt = odefun( t , y )
    Tr = ctte;
    L = ctte;
    dydt = sin(y) + t - y/t - Tr + L*L/Tr; % <<=== à modifier en fonction de ta première étape
end

... et tu l'enregistres sous le fichier odefun.m

Troisième étape : tu as enfin le droit d'écrire :

Code :

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

[T,Y] = ode45(@odefun,tspan,y0)

(Attention à ne pas oublier le @)

Si tu butes, réponds juste à la première étape, stp.

Cordialement,

[cpalperou]

Bonjour et merci à vous: j'ai réussi!!

voici ce que j'ai fait:

Code :

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

1
2
3
function dy=rigid(t,y,f,ft)
f=interp1((ft,f,t)
dy = f;

puis dans la fenêtre principale de MATLAB:

Code :

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

1
2
3
4
5
f=1./(L.*Tr)-1./L;
Tspan=20:120;
ft=Tspan;
Ic=1e-2;
[T Y] = ode45(@(t,y) rigid(t,y,f,ft),Tspan,Ic)

et bingo!
Par contre, Tachmou, il me semble que l'exemple 3 de la documentation (lien fourni par DUT) te contredit quand tu affirmes "qu'il est impossible de mettre autre chose que t et y en entrée de cette fonction" !

[Invité]

Code :

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

1
2
f=interp1((ft,f,t)
dy = f;

ou directement dy = interp1((ft,f,t)

Par contre, Tachmou, il me semble que l'exemple 3 de la documentation (lien fourni par DUT) te contredit quand tu affirmes "qu'il est impossible de mettre autre chose que t et y en entrée de cette fonction" !

Ce n'est pas tout à fait vrai, la syntaxe @(t,y) rigid(t,y,f,ft) crée une fonction anonyme équivalent à écrire ;

Code :

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

1
2
function r = fun(t,y)
r = rigid(t,y,f,ft);

Avec l'avantage de prendre les variables f et ft dans le workspace courant.
Tu ne fournis au final qu'une fonction ne prenant que deux arguments t et y, et c'est bien la seule forme qu'acceptera les fonctions ode.

[tachmou]

+1 avec Winjerome

C'est pas simple à comprendre je le conçois... mais ça fait plaisir de voir que tu remets en question : ça montre que tu te pose les bonnes questions ! Donc +1 à toi aussi

A+