Bonjour la compagnie,

Voilà, j'aimerais faire le tp suivant : http://www.thphys.may.ie/CompPhysics...pingFaucet.pdf

Voilà, J'aimerais afficher le diagramme des bifurcations avec Matlab et j'ai écris un programme mais il ne marche pas comme il faudrait:

Code:
clc;clf;
% constantes
rho = 1; % masse volumique de l'eau [g/cm^3]
xc = 0.25; % hauteur critique [cm]
alpha = 0.5; % coeff proportionnalite [s/m]
% conditions initiales
x0 = 0; v0 = 0; m0 = 0.01;
for debit = 0.57 : 0.001 : 0.7 % fenetre du diagramme de bifurcations
    debit % pour voir où en est le programme
    Niter = 200; % nombre d'iterations
    goutte = zeros(Niter,4); % [tc xc vc mc]
    goutte(1,:) = [0 x0 v0 m0];
    for k = 2 : Niter
        [t,x] = ode45('robinet',[0 2],[x0 v0 m0 debit]);
        % recherche point critique de detachement
        rc = min(find(x(:,1)>xc)); % Rang Critique
        tc = t(rc,1);
        Xc = x(rc,1); % hauteur critique
        Vc = x(rc,2); % vitesse critique
        Mc = x(rc,3); % masse critique
        % detachement de la goutte
        deltaM = alpha * Mc * Vc; % masse de la goutte tombante
        rgt = (3*deltaM/(4*pi*rho))^(1/3); % Rayon de la Goutte Tombante
        % nouvelles valeurs
        x0 = Xc - (rgt * deltaM / Mc);
        m0 = Mc - deltaM;
        v0 = Vc;
        goutte(k+1,:) = [tc x0 v0 m0];
    end;
    % dessin des 50 dernieres iterations
    subplot(1,2,1);plot(ones(50,1)*debit,goutte(end-49:end,1),'b.','markersize',1); hold on;
    subplot(1,2,2);plot(ones(50,1)*debit,goutte(end-49:end,3),'b.','markersize',1); hold on;
end;
avec
Code:
function dx = robinet(t,x)
% constantes 
g = 980; % cste de gravite [cm/s^2]
k = 475; % raideur du ressort [dyn/cm]
b = 1; % coeff de frottement visqueux [g/s]
dx = zeros(4,1); % [position; vitesse; masse]
% equation
dx(1) = x(2); % dx/dt = v;
dx(2) = g - (x(1) * k/x(3))  - x(2)*(b+x(4))/x(3); % dv/dt = g - (k/m)*x - v*(b+Q)/m;
dx(3) = x(4) ; % dm/dt = Q;
Voilà, ça affiche pas ce je voudrait et je ne sais pas pourquoi. Alors si quelqu'un sait comment faire, ce serait cool