// TP4 Correction : Résolution d'équations non linéaires, Introduction aux EDOs clear // Exercice 1 function y=f(x) y=x-1/5*sin(x)-1/2; endfunction //question 1; On a besoin de la fonction f pour le critère d'incrément et de la fonction phi pour l'itération function y=phi(x) y=1/5*sin(x)+1/2 endfunction function [k,x]=pointfixe(alpha,tol,g,h) x=alpha;k=0 while norm(h(x))>tol x=g(x); k=k+1; end endfunction //test : alpha=-5; tol=1e-5; [kfix1,xfix1]=pointfixe(alpha,tol,phi,f) tol=1e-9; [kfix2,xfix2]=pointfixe(alpha,tol,phi,f) tol=1e-12; [kfix3,xfix3]=pointfixe(alpha,tol,phi,f) //question 2 : on a besoin pour la mise en place de la méthode de Newton de la fonction f que l'on a déjà et de sa dérivée. On la calcule à la main et on l'implémente. function y=df(x) y=1-1/5*cos(x); endfunction function [k,x]=newton(alpha,tol,g,dg) x=alpha; k=0; while norm(g(x))>tol x=x-g(x)./dg(x); k=k+1; end endfunction alpha=-5;tol=1e-5; [knew1,xnew1]=newton(alpha,tol,f,df) tol=1e-9; [knew2,xnew2]=newton(alpha,tol,f,df) tol=1e-12; [knew3,xnew3]=newton(alpha,tol,f,df) //La méthodede Newton converge beaucoup plus vite : plus la tolérance est petite, plus la différence entre Newton et la méthode de points fixe est flagrante //question 3 //Pour cette question on modifie les fonctions ci dessus pour éviter le test de tolérance et pouvoir calculer le nombre d'éléments voulu de la suite. function [err,resx]=pointfixe2(alpha,K,g,h) err=[norm(h(alpha))]; resx=[alpha]; x=alpha; for k=1:K x=g(x); resx=[resx,x]; err=[err norm(h(x))]; end endfunction function [err,resx]=newton2(alpha,K,g,dg) x=alpha;resx=[alpha]; err=[abs(g(alpha))]; for k=1:K x=x-g(x)./dg(x); resx=[resx x]; err=[err abs(g(x))]; end endfunction alpha=-5 K=15; [err1,resx1]=pointfixe2(alpha,K,phi,f) [err2,resx2]=newton2(alpha,K,f,df) plot(0:K,log(err1),0:K,log(err2)) legend('ptfixe','newton') //---------------------------------------------------------------------- //Exercice 2 : Un peu de géométrie //question 1 : Les solutions sont (sqrt(3/2),sqrt(1/2)), (sqrt(3/2),-sqrt(1/2)), (-sqrt(3/2),sqrt(1/2)), (-sqrt(3/2),-sqrt(1/2)) //A noter que l'on peut tracer ces courbes comme suit : x1=linspace(-sqrt(2),sqrt(2),50); y1=sqrt(2-x1.^2); z1=-sqrt(2-x1.^2); x2=linspace(-10,-1,500); y2=sqrt(x2.^2-1); z2=-sqrt(x2.^2-1); x3=linspace(1,10,500); y3=sqrt(x3.^2-1); z3=-sqrt(x3.^2-1); scf()//cette commande permet d'ouvrir une nouvelle fenetre de graphiques plot(x1,y1,x1,z1,x2,y2,x2,z2,x3,y3,x3,z3) plot(,sqrt(3/2),sqrt(1/2),'o',-sqrt(3/2),sqrt(1/2),'o',sqrt(3/2),-sqrt(1/2),'o',-sqrt(3/2),-sqrt(1/2),'o') //question 2 function Y=systeme(X) x=X(1); y=X(2); Y=[x^2+y^2-2;x^2-y^2-1]; endfunction function A=Diffsysteme(X) x=X(1);y=X(2); A=[2*x 2*y;2*x -2*y]; endfunction function [X,comp]=Newton2D(X0,tol) X=X0;Xaux=X0;err=1;comp=0 while err>tol Y=systeme(X); A=Diffsysteme(X); X=X-inv(A)*Y; err=norm(X-Xaux); Xaux=X; comp=comp+1; end endfunction X0=[1;-1]; tol=1e-9; [X1,comp1]=Newton2D(X0,tol); X0=[1;1]; tol=1e-9; [X2,comp2]=Newton2D(X0,tol); X0=[-1;1]; tol=1e-9; [X3,comp3]=Newton2D(X0,tol); X0=[-1;-1]; tol=1e-9; [X4,comp4]=Newton2D(X0,tol); // On pourrait aussi s'intéresser à la suite des points obtenus par la méthode de Newton pour les voir se rapprocher des solution, ou à l'évolution de l'erreur d'incrément en fonction du nombre d'itération... function [resX,X,err]=Newton2D(X0,K) X=X0;resX=[X0];Xaux=X0,err=[]; for k=1:K Y=systeme(X); A=Diffsysteme(X); X=X-inv(A)*Y; resX=[resX, X]; err=[err,norm(X-Xaux)]; Xaux=X; end endfunction K=15;X0=[1;-1]; [resX,X,err]=Newton2D(X0,K) x1=linspace(0,sqrt(2),50); z1=-sqrt(2-x1.^2); x3=linspace(1,2,500); z3=-sqrt(x3.^2-1); scf() plot(x1,z1,x3,z3) plot(sqrt(3/2),-sqrt(1/2)) plot(resX(1,:),resX(2,:),'*') scf() plot(1:15, log(err)) //------------------------------------------------------------------- //Exercice 3 : Pendule simple, introduction aux EDOs // //question 1 : la solution explicite est donnée par // theta(t)=theta0*cos(t)+theta1*sin(t)=> x(t)=theta0 cos(t)+theta1 sin(t) //et y(t)=-theta0 sin(t)+theta1 cos(t). // //question 2 function Y=pendule(X)//fonction qui permet de décrire l'EDO en vectoriel x=X(1);y=X(2); Y=[y;-x]; endfunction function y=exact(theta0,theta1,t)//solution exacte y1=theta0*cos(t)+theta1*sin(t); y2=-theta0*sin(t)+theta1*cos(t); y=[y1;y2]; endfunction function [discrt,resY]=Eulerex(Y0,h,T)//méthode d'Euler Explicite vectorielle //discrt est la discrétisation en temps //resY est la liste des valeurs de y_k sur cette discrétisation N=floor(T/h); discrt=[0];t=0; resY=[Y0]; Y=Y0; for k=1:N Y=Y+h*pendule(Y); resY=[resY, Y]; t=t+h; discrt=[discrt,t]; end endfunction function [discrt,resY]=Heun(Y0,h,T)//comme pour Euler avec mais pour Heun N=floor(T/h); discrt=[0];t=0; resY=[Y0]; Y=Y0; for k=1:N Y1=Y+h*pendule(Y) Y=Y+h/2*(pendule(Y)+pendule(Y1)); resY=[resY, Y]; t=t+h; discrt=[discrt,t]; end endfunction theta0=0.01; theta1=0; Y0=[theta0;theta1]; h=1e-3; T=50; [t,resYE]=Eulerex(Y0,h,T) [discrt,resYH]=Heun(Y0,h,T) Z=exact(theta0,theta1,t) scf()//tracé des angles au cours du temps : première ligne des vecteurs obtenus plot(t,resYE(1,:),t, resYH(1,:),t,Z(1,:)) legend('Euler ex','Heun','exact') title('Angle en fonction du temps') scf()//tracé des vitesses angulaires au cours du temps : 2eme ligne des vecteurs obtenus plot(t,resYE(2,:),t,resYH(2,:),t,Z(2,:)) legend('eulerex','Heun','exact') title('vitesse angulaire en fonction du temps') scf() //plan de phase plot(resYE(1,:),resYE(2,:),resYH(1,:),resYH(2,:)) legend('EulerE','Heun') plot(theta0,theta1,'o') title('plan de phase') scf() //energie conservée ? energieEuler=resYE(1,:).^2+resYE(2,:)^2; energieHeun=resYH(1,:).^2+resYH(2,:)^2; energieExact=Z(1,:)^2+Z(2,:)^2; plot(t,energieEuler,t,energieHeun,t,energieExact) legend('euler','Heun','exact') title('energie au cours tu temps')