// TP 3 MT12 P25 /// /// /// /// /// /// /// /// PARTIE 1 Intégration d'une fonction continue sur un intervalle borne clear // ------------ // Q1 //------------ function [Imoins, Iplus, N] = Riemann(f, a, b, epsilon) nsubdiv = 20 // Voir la Remarque juste après N = 1 done = %F while (~done) N = 2 * N Imoins = 0 Iplus = 0 for k = 1:N Ik = [ a+(k-1)*(b-a)/N, a+k*(b-a)/N ] // (1) subdiv = linspace(Ik(1), Ik(2), nsubdiv) // Subdivision de l’intervalle fmin = min(f(subdiv)) fmax = max(f(subdiv)) Imoins = Imoins + fmin Iplus = Iplus + fmax end Imoins = Imoins * (b-a)/N Iplus = Iplus * (b-a)/N done = (Iplus - Imoins) < epsilon end // (while) endfunction function y=sinus_pi(x) y=sin(%pi.*x) endfunction [Imoins, Iplus,N] = Riemann(sinus_pi, 0, 1, 10^-4); // Integrale int_0^1 sin(pi t) de = -cos( pi t )/pi ]_0^1= 2/pi z_sin= integrate('sinus_pi(x)', 'x', 0, 1); disp("fonction TP: ", Imoins, Iplus) disp("fonction scilab integrate: ",z_sin) disp("valeur exacte (calculée): ", 2/%pi) // ------------ // Q2 //------------ // Plot FONCTION sin(pi x) sur 0, 1 // Intégrale de Riemann zx=linspace(0,1,200) zy=sinus_pi(zx) scf() plot(zx,zy) title("fonction x ---> sin (pi x)") // Méthode des rectangles à gauche avec borne sup // ------------ // Q3 //------------ // Avec Borne sup function y=der_sin_pi(x) // fonction dérivée y=%pi.*cos(%pi.*x) endfunction // valeur de l'intégrale, son approximation, leur différence et la borne sup function [I,I_N, dif_I, bor_sup] = Rectangle_g(f, a, b,N) nsubdiv = 20; // deb: borne sup x_sup=linspace(a,b,nsubdiv) y_sup=abs(der_sin_pi(x_sup)) bor_sup=max(y_sup).*((b-a).^2./(2*N)) // fin : borne sup x=linspace(a,b,N) // deb : methode rectangle I=integrate('sinus_pi(x)','x',a,b) I_N = 0; for k = 1:N xk = a+(k-1)*(b-a)/N; I_N = I_N + sinus_pi(xk); end I_N = I_N.*((b-a)./N); dif_I=abs(I_N-I); endfunction [I,I_N, dif_I, bor_sup] = Rectangle_g(sinus_pi, 0, 1,100) // Sans La borne sup function [I_rec, I_diff, N] = Rectangle_alt(f, a, b, epsilon) N = 1 done = %F while (~done) N = 2 * N I_rec = 0 for k = 1:N xk= a+(k-1)*(b-a)/N I_rec = I_rec + f(xk) end I_rec= I_rec * (b-a)/N I=integrate('f(x)','x',a,b) I_diff=I - I_rec done = I_diff< epsilon end // (while) endfunction [I_rec, I_diff, N] = Rectangle_alt(sinus_pi, 0, 1, 10^(-4)) //////////////// PARTIE 2 CONVERGENCE PONCTUELLE SERIE DE FOURIER FONCTION PERIODIQUE REGULIERE // ------------ // Q4 //------------ //Graphe de la fonction sur plusieurs périodes de longueur 2 function y=Perio_cont(x) y= (abs(x)-(2.*int(abs(x)/2))).*(modulo(floor(abs(x)),2)==0) + (2-(abs(x)-(2.*int(abs(x)/2)))).*(modulo(floor(abs(x)),2)==1) endfunction // Explication //1. fonction paire donc f(-x)=f(x) donc resultat en fonction de abs(x) ; //2. quand x in [ 2k, 2k+1[, k in N alors f(x) = y pour y = x- 2* int(x/2) avec int(x/2)= quotient de la div Eucli de x par 2 ; //condition sur x : ""reste de la div Eu de la partie entiere de x par 2 = 0" et f(x) = x- 2*(quotient de la div euclidienne de x par 2) // 3 quand x in [ 2k+1 , 2(k+1) ] k in N alors f(x) = 2- y (pour se ramener à une valeur entre 0 et 1) // condition sur x : "(modulo(floor(abs(x)),2)==1" scf() zx=linspace(-6,6,601) zy=Perio_cont(zx); title('fonction triangulaire') plot(zx,zy) // ALT : Graphe de la fonction uniquement sur [0,2] function y=Perio_cont_alt(x) y= x.*(0<=x & x<1 )+(2-x).*(x>=1 & x<=2) endfunction scf() zx=linspace(0,2,201) zy=Perio_cont_alt(zx) title('fonction triangulaire') plot(zx,zy) // ------------ // Q5 //------------ // Coefficients de Fourier à la main // a0 =1 // bk = 0 (fonction paire) // ak = -(1-(1)^k)*(2/(k pi)^2) pour k >0 // et donc ak = 0 si k>0 pair // ak = -4/((k pi)^2) si k impair // on définit fN = a0/2 + sum_k>=1^N (ak cos( k pi x) + bk sin(k pi x)) // ------------ // Q6 //------------ // On simplifie alors f_(2N+1)(x)) par 1/2 - 4/(pi^2)sum_{l=0}^N (cos(pi (2l+1) x) )/((2l+1)^2) // ------------ // Q7 //------------ function y=somme_partielle_Fourier_paire(N,x) dL=0 for l=0:N k=2.*l+1; // le k impair du ak non nul dL=dL+ cos(%pi.*k.*x)/(k.^2); end y=(1/2)-(4/%pi^2).*dL; endfunction zx=linspace(0,2,201) ; zy=Perio_cont(zx); N=[2, 4, 8, 16, 32, 64, 128,256]; zN=zeros(length(N),length(zx)); for i=1:length(N) zN(i,:)=somme_partielle_Fourier_paire(N(i),zx); end zz=[zy', zN']; zx=zx'; scf() plot(zx,zz,LineSpec=1:9) legend(['True','N=2','N=4','N=8','N=16','N=32','N=64','N=128','N=256']) ////////////////// PARTIE 3. "Phénomène de GIBBS POUR UNE FONCTION PERIODIQUE DISCONTINUE " function y=Perio_Discont_supp(x) // uniquement sur le support [-1,1] y= (-1-x).*(-1 <=x & x <0) + (1-x).*(0<=x & x<= 1) endfunction zx=linspace(-1,1,201) ; zy=Perio_Discont_supp(zx) // ------------ // Q8 //------------ scf() plot(zx,zy) // ------------ // Q9-Q10 //------------ // a_n(f)=0 fonction impaire // b_n(f)=1/(pi n) function y=somme_partielle_Fourier_Impaire(N, x) dL=0 for k=1:N dL=dL+2/(k*%pi).*sin(k*%pi.*x) end y=dL endfunction zx=linspace(-1,1,101) ; zy=Perio_Discont_supp(zx); N=[2, 4, 8, 16, 32, 64, 128,256] zN=zeros(length(N),length(zx)); for i=1:length(N) zN(i,:)=somme_partielle_Fourier_Impaire(N(i),zx); end zz=[zy', zN']; zx=zx'; scf() plot(zx,zz,LineSpec=1:9) legend(['True','N=2','N=4','N=8','N=16','N=32','N=64','N=128','N=256'],4) // lower right // ------------ // Q11 //------------ // Value of fN(0)=1 forall N et f(0)=1 "fN (0) - f(0) " zy0N=zeros(length(N),1); dif0N=zeros(length(N),1); for i=1:length(N) zyN0(i)=somme_partielle_Fourier_Impaire(N(i),0); dif0N(i)=zyN0(i)-Perio_Discont_supp(0); end dif0N // ------------ // Q12 //------------ // PAS DE CONVERGENCE PONCTUELLE EN x=0 // c'est cohérent avec le cours // ------------ // Q13 //------------ // Illustration de la cv dans L2 ; on approche la norme L2 au carré de la différence function y=JKN(K,N) y=0 for k=0:K-1 y=y+(abs(Perio_Discont_supp( -1+k*(2/K)) - somme_partielle_Fourier_Impaire(N,-1+k*(2/K))))^2 end y=(2/K)*y; endfunction N=[ 64, 128, 256, 512, 1024, 2048 ] JKN_vect=zeros(length(N),1); K=600 for i=1:length(N) JKN_vect(i,:)=JKN(K,N(i)); end JKN_vect // JKN tend vers 0 avec N qui grandit ce qui est cohérent avec le cours .