Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

  1. #1
    Candidat au Club
    Problème Résolution système équation différentiel en utilisant ODE45
    Bonsoir à tous,
    je souhaiterais résoudre sur Matlab ce système d'équation :



    J'ai essayé le code ci dessous :

    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
     
    function ypoint = ab(t,y) 
    global mumax1 mumax2 ks1 ks2 ki2 k1 k2 k3 k4 
    ypoint(1) = (mumax1*y(3)*y(1)/(y(3)+ ks1)); 
    ypoint(3) = -k1*(mumax1*y(3)*y(1)/(y(3)+ks1));
    ypoint(2) = (mumax2*y(4)*y(2)/(y(4)+ks2+(y(4)^2/ki2)))
    ypoint(4) = k2*(mummax1*y(3)*y(1)/(y(3)+ks1))-k3*(mumax2*y(4)*y(2)/(y(4)+ks2+(y(4)^2/ki2)));
    ypoint(5)  = k4*(mumax2*y(4)*Y(2)/(y(4)+ks2+(y(4)^2/ki2)));
    ypoint = ypoint(<img src="images/smilies/icon_smile.gif" border="0" alt="" title=":)" class="inlineimg" />;
    end
     
    global mumax1 mumax2 ks1 ks2 ki2 k1 k2 k3 k4 q
     mumax1=0.4; mumax2= 0.4; ks1= 35; ks2=4; ki2= 170; k1= 50; k2= 50; k3= 15 k4= 75;% paramètres
    tfinal= 100; % temps final
    y01= 0.4; y02= 0.01; y03= 10; y04= 2; y05= 0;% conditions initiales
    [t,y]= ode45('ab',[0 tfinal],[y01 y02 y03 y04 y05])% résolution
    y1 = y(:,1);%extractions de Y1 Y2 Y3 Y4 et y5 
    y2 = y(:,2);
    y3 = y(:,3);
    y4 = y(:,4);
    y5 = y(:,5);
    q = k4*(mummax2*y4*y2/(y4+ks2+(y4^2/ki2)));
    plot(t,y1,'r') % y1 fonction de temps
    plot(t,y2,'g') % y2 fonction de temps
    plot(t,y3,'o') % y3 fonction de temps
    plot(t,y4,'a') % y4 fonction de temps
    plot(t,y1,t,y2,t,y3,t,y4)% y4 fonvtion de temps
    plot(t,y5)% y4 fonvtion de temps
    plot(t,q,'z');
    plot(t,y5,'r',t,q,'p');


    Mais je recois cette erreur :

    >>Untitled2
    Error : File : Untitled2.m Line : 10 column : 1
    This statement is not inside any function.
    (It follows the END that terminates the definition of the function "ab".)


    Pouvez vous m'aider svp

  2. #2
    Membre actif
    Bonjour,

    La définition de la fonction ab doit se trouver à la fin du script, et non au début.
    Si tu n'es pas familier de la fonction ode45, je te conseille de commencer par tenter de résoudre une EDO plus simple, voire de te baser sur un example publié par Mathworks. e.g. https://fr.mathworks.com/help/matlab...ref/ode45.html
    Note par ailleurs qu'il serait préférable d'éviter les variables globales. Mathworks propose diverses manière de procéder : https://fr.mathworks.com/help/matlab...5.html#bu3uhuk

  3. #3
    Futur Membre du Club
    Citation Envoyé par Antoane Voir le message
    Bonjour,

    La définition de la fonction ab doit se trouver à la fin du script, et non au début.
    Si tu n'es pas familier de la fonction ode45, je te conseille de commencer par tenter de résoudre une EDO plus simple, voire de te baser sur un example publié par Mathworks. e.g. https://fr.mathworks.com/help/matlab...ref/ode45.html
    Note par ailleurs qu'il serait préférable d'éviter les variables globales. Mathworks propose diverses manière de procéder : https://fr.mathworks.com/help/matlab...5.html#bu3uhuk
    bonsoir,
    je voudrais que vous m'aider; résoudre cette équation différentielle à l'aide de runde kutta 4
    voici le code de mon programme de la fonction de l'autre coté de l'égalité
    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
    function [g]= foc(t,y)
    gama=10; f=1500; alpha=0.01; Mo=14; V=1.4e-6;
    tg=alpha/gama;te=alpha/Mo*V^2; b=10e2; d=10e8;
    w=2*3.14*f;
    A=[0.5-3*d*V^2/16 0 -0.25+d*V^2/8 0 -d*V^2/32];
    B=[0 -0.25*b*V 0 b*V/12 0];
    % t= linspace(0,20,1000) ;
    s=0;
    e=V*sin(w*t);
    W=Mo*(0.5-e*b/3-d/4*e.^2).*e.^2;
    figure(1)
    plot(t,W)
    for n=2:5;
        s=s+A(1,n)*cos(n*w.*t)+B(1,n).*sin(n*w.*t);
    end
    r=0.5*A(1,1)+s;
    % g=W/(te*Mo*V^2)-y/tg;
    g=1/te*r-y/tg;
    end

    voici mon programme runge kutta 4
    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
    function [y,t] = RK4(f,tmin,tmax,Nint,y0) % Méthode de Runge Kutta d’ordre 4
    % Nint - nombre de sous intervalles
    % tmin - temps t0
    % tmax - temps t0 + T
    % f est une fonction avec comme arguments t et y : f(t,y(t))
    % y0 contient les valeurs des conditions limites
    h = (tmax-tmin)/Nint ; % valeur du pas (G)
    t = linspace(tmin,tmax,Nint+1) ; % vecteur de t discrétisé t=[tmin,tmax]
    n = length(t);
    y = zeros(1,n); % on fixe la taille du tableau y
    y(1,1) = y0 ;
    t(1,1)=0;% y contient les solutions de y(tn)n = 1, ..., Nint + 1
    for n = 2:Nint-1
     t(n+1)=t(n)+h;
     k1=h*f(t(n),y(n));
     k2=h*f(t(n)+(h/2),y(n)+(k1/2));
     k3=h*f(t(n)+(h/2),y(n)+(k2/2));
     k4=h*f(t(n)+h,y(n)+k3);
    y(n+1)=y(n)+(k1/6)+(k2/3)+(k3/3)+(k4/6);
    end 
    end


  4. #4
    Membre actif
    Bonjour,

    Et alors ?
    Comment ce code est-il exécuté ? Qu'est ce que ça donne ? Quel éventuel message d'erreur ?
    etc.

###raw>template_hook.ano_emploi###