/* TP 6 Résolution approchée d'équations différentielles */ /* Partie 2 Résolution numérqiue d'EDO */ /* Q1 : solution y(t)=10 exp(-50t) /* /* Q2 solutions approchées */ function y=sol(t,Y0) y=Y0.* exp(-50.*t) endfunction function y=f(t) y= - 50.*t endfunction function j=I(i) j = i+1; endfunction /* Solution approchée par Euler */ function Y=EulerExpl(Y0,T,N,f) dt=T/N; Y=zeros(length(Y0),N+1); Y(:,I(0)) =Y0 ; for i=0:(N-1) Y(:,I(i+1))=Y(:,I(i)) + dt.*f(Y(:,I(i))) end endfunction /* Solution approchée par Heun */ function Y=Heun(Y0,T,N,f) dt=T/N; Y=zeros(length(Y0),N+1); Y(:,I(0)) =Y0 ; for i=0:(N-1) Y(:,I(i+1))=Y(:,I(i)) + (dt/2).*( f(Y(:,I(i))) + f(Y(:,I(i)) + dt.*f(Y(:,I(i))))) end endfunction /* Solution approchée par RK4 */ function Y=RK4(Y0,T,N,f) dt=T/N; Y=zeros(length(Y0),N+1); Y(:,I(0)) =Y0 ; for i=0:N-1 k1=f(Y(:,I(i))) k2=f( Y(:,I(i)) + k1.*(dt/2)) k3=f(Y(:,I(i))+ k2.*(dt/2)) k4=f(Y(:,I(i)) + k3.*dt) Y(:,I(i+1))=Y(:, I(i)) + (dt/6).*(k1+2.* k2+2*k3+k4); end endfunction /*Application */ T=1 ; N=50; Y0=10; dt=T/N; t=0:dt:T; ysol=sol(t,Y0); yappE=EulerExpl(Y0,T,N,f); yappH=Heun(Y0,T,N,f); yappRK4=RK4(Y0,T,N,f); scf(1) subplot(3,1,1) plot(t,ysol) plot(t,yappE,'g+-') title("Methode Euler") legend(["sol exacte", "sol approchée"]) subplot(3,1,2) plot(t,ysol) plot(t,yappH,'g+-') title("Methode Heun") legend(["sol exacte", "sol approchée"]) subplot(3,1,3) plot(t,ysol) plot(t,yappRK4,'g+-') title("Methode Runge-Kutta ordre 4") legend(["sol exacte", "sol approchée"]) /*Q3 Réponse : la troisième méthode */ /* Systeme de 2 EDO découplées*/ /* Q4 En posant y(t)=(y1(t),y2(t)) , f (y1(t),y2(t)) = A () =(y1(t),y2(t)), avec la matrice carré 2*2 A= (-50t, 0 ; 0, 3t)) y(0)=(10,2) Solution exacte y1(t)=10 exp(-50t) y2(t)=2 exp(3t) */ /*Q5 */ function y=solD(t,Y0) y=[Y0(1).*exp(-50.*t) ; Y0(2).*exp(3.*t)] endfunction function y=fD(Y) y= [- 50.*Y(1,:);3.*Y(2,:)] endfunction T=1 N=50 Y0=[10,2] TT=1/N; t=0:TT:1 ysolD=solD(t,Y0); yappE=EulerExpl(Y0,T,N,fD); yappH=Heun(Y0,T,N,fD); yappR=RK4(Y0,T,N,fD) ; /* Graphe des solutions du système d'EDO découplées */ scf(2) subplot(3,1,1) plot(t,ysolD(1,:)) plot(t,yappE(1,:),'g+-') plot(t,ysolD(2,:)) plot(t,yappE(2,:),'r+-') title("Methode Euler") legend(["sol exacte (1)", "sol approchée (1)", "sol exacte (2)", "sol approchée (2)"],[2]) subplot(3,1,2) plot(t,ysolD(1,:)) plot(t,yappH(1,:),'g+-') plot(t,ysolD(2,:)) plot(t,yappH(2,:),'r+-') title("Methode Heun") legend(["sol exacte (1)", "sol approchée (1)", "sol exacte (2)", "sol approchée (2)"],[2]) subplot(3,1,3) plot(t,ysolD(1,:)) plot(t,yappR(1,:),'g+-') plot(t,ysolD(2,:)) plot(t,yappR(2,:),'r+-') title("Methode Runge-Kutta ordre 4") legend(["sol exacte (1)", "sol approchée (1)", "sol exacte (2)", "sol approchée (2)"],[2]) /* Partie 2.3 Modèle Proies-Prédateurs*/ /* Q6*/ function ans=fPP(X) ans=[3*X(1,:)- X(1,:).*X(2,:) ;-2*X(2,:) + X(1,:).*X(2,:)] endfunction /* Q7-- Q8*/ N=5000 T=20 Y0=[1,2] TT=1:N t=(TT.*T)/N; t=[0,t]; yappEPP=EulerExpl(Y0,T,N,fPP); yappHPP=Heun(Y0,T,N,fPP); yappRPP=RK4(Y0,T,N,fPP); scf(4) subplot(3,1,1) comet(yappEPP(1,:),yappEPP(2,:)) title("Solution approchée Methode Euler") subplot(3,1,2) comet(yappHPP(1,:),yappHPP(2,:)) title("Solution approchée Methode Heun") subplot(3,1,3) comet(yappRPP(1,:),yappRPP(2,:)) title("Solution approchée Methode Runge-Kutta ordre 4") /* Partie 3 Résolution numérique de l’équation de la chaleur sur domaine non borné */ /* Q9 -- Q10 */ // discretisation domaine spatial L=100; // on considere le carre [-50,50] M=100; dx=L/M;// pas d'espace x=-L/2:dx:L/2-dx; // discretisation domaine frequentiel xi= (1/L)*[-M/2:M/2-1]'; xi=fftshift(xi);// utile pour la suite // discretisation domaine temporel T=20; dt=1e-3;// pas de temps t=0:dt:T; N=length(t); /* Fonction f v: M lignes et 1 colonne */ function res=fCh(v,xi,a) res=-a.*4.*%pi^2.*xi.^2.*v; endfunction a=1 /* ATTENTION : pour utiliser RK4, il est nécessaire de déclarer "f" comme fonction d''un seul argument cela pré-suppose de définir à l'extérieiur de la fonction le xi et le a ... pas très propre d'un point de vue implémentation d'une fonction ... */ function res=fChaleur(v) res=fCh(v,xi,a); endfunction /* u0 */ u0=zeros(1,length(x)); u0((L/2-L/10)/dx:(L/2+L/10)/dx)= 1; // Indice sur les x qui satisfont la condition u0=u0'; /* Q11 */ a=1; Y0=fft(u0); uhat=RK4(Y0,T,N,fChaleur); /* Q12*/ for k=1:N u(:,k)=ifft(uhat(:,k)); end clf() scf(5) plot(x,u(:,1)',"black") plot(x,u(:,floor(N/10)),"g") plot(x,u(:,floor(N/2))',"r") plot(x,u(:,N),"b") legend("t=0","t=2","t=10","t=20")