Résoudre une EDO 2ème orde avec MATLAB

Résoudre une EDO 2ème orde avec MATLAB - Divers - Programmation

Marsh Posté le 02-03-2005 à 16:39:13    

bon voila ca marche pas ?????
c'est un système de deux equations du 1er orde :
 
dV/dt = W
dW/dt = -(L+RC^2)*W/(LCR)-2*V/(LC)+e/(LC)
W(0) = V(0) = 0
 
La méthode que je dois utiliser est RK4, ci-dessous. Mais ca marche pas, le graphique obtenut est faut.
 
 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
%*************************************
% FICHIER DE COMMANDE "probleme2.m"
%*************************************
 
clear all;  
 
%PARAMÈTRES NUMÉRIQUES
%Inductance [H]
L = 0.1;
%Capacité [F]
C = 0.001;
%Résistance [ohm]
R = 10;
%Tension d'entrée [V]
e = 5;
 
%INITIALISATIONS
%thermes pour l'equation Y2
a = -(L+R*C^2)/(L*C*R);
b = -2/(L*C);
c = e/(L*C);
%indice de temps
i = 1;
%Vecteur temps
t(1) = 0;
%Pas de discrétisation
h = 0.01;
%Conditions initiales
y(:,1) = [0;0];
 
%SOLUTION APPROCHEE À PAS DE DISCRETISATION CONSTANT
%=======================================================================
%calculée par la méthode de RUNGE-KUTTA d'ordre 4
%-----------------------------------------------------------------------
 
tic;
 
test = 1;
while i < 50 %test > 0.0001 %boucle jusqu'au regime permanant
   
   % f(wn,tn)
   k1 = fct(y(:,i),t(i),a,b,c);
   
   % f(wn+h*k1/2,tn+h/2)  
   k2 = fct(y(:,i)+h*k1/2,t(i)+h/2,a,b,c);
   
   % f(wn+h*k2/2,tn+h/2
   k3 = fct(y(:,i)+h*k2/2,t(i)+h/2,a,b,c);
   
   % f(wn+h*k3,tn+1)
   k4 = fct(y(:,i)+h*k3,t(i)+h,a,b,c);
   
   y(:,i+1) = y(:,i)+(k1+2*k2+2*k3+k4)/6;
   t(i+1) = t(i) + h;
   i = i+1;
   
   %test = y(1,i)-y(1,i-1);
     
end
toc  
 
 
%AFFICHAGE GRAPHIQUE 2D
%=======================================================================
%Trace de v(t) la tension du circuit  
%-----------------------------------------------------------------------
figure(2)
plot(t,y(1,:),'r','LineWidth',2);
grid on;
xlabel('Temps (en s)');
ylabel('Tension (en V)');  
title('Evolution de v(t)');
%axis([0 max(t) 0 max(v)+5])
%text(max(t)/2,400,['Fusion du fusible à ',num2str(max(t)),' secondes'],'BackgroundColor',[.7 .9 .7])
 
figure(3)
plot(t,y(2,:),'LineWidth',2);
grid on;
xlabel('Temps (en s)');
ylabel('Tension (en V)');  
title('Evolution de w(t)');
 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
%******************************************************                      
%FICHIER de FONCTION "fct.m"                                              
%******************************************************                      
                                                                             
                                                   
function f = fct(y,t,a,b,c);                                                          
 
f = [y(2);a*y(2)+b*y(1)+c];
 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
 
Bon pour ne pas que vous aillez trop de truc à faire pour debug mon programme ci-dessous il y a le système résolu avec ODE45.
 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
 
clear all;
 
tic;
[T,Y] = ode45(@fct_ode,[0 0.25],[0 0]);
toc
 
figure(1)
plot(T,Y(:,1),'r','LineWidth',2);
%axis([0 120 -1.5 1.5])
grid on;
 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
function dy = fct(t,y)
 
%PARAMÈTRES NUMÉRIQUES
%Inductance [H]
L = 0.1;
%Capacité [F]
C = 0.001;
%Résistance [ohm]
R = 10;
%Tension d'entrée [V]
e = 5;
 
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = -(L+R*C^2)*y(2)/(L*C*R)-2*y(1)/(L*C)+e/(L*C);
 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
 
Voila merci de m'aider.
 
Pierre.

Reply

Marsh Posté le 02-03-2005 à 16:39:13   

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed