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.