From: couchot Date: Mon, 9 Sep 2013 07:50:33 +0000 (+0200) Subject: init X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/cours-mesi.git/commitdiff_plain/fe25b30b852e5ee9b37ad770d19f3b1a070349e3?ds=inline init --- diff --git a/partiels/main.tex b/partiels/main.tex new file mode 100644 index 0000000..57c391d --- /dev/null +++ b/partiels/main.tex @@ -0,0 +1,117 @@ +\documentclass[10pt,a4paper,french]{article} +\usepackage[francais]{babel} +\usepackage[utf8]{inputenc} +\usepackage{a4} +\usepackage{amsmath} +\usepackage{amsfonts} +\usepackage{amssymb} +\usepackage{framed} +\usepackage{dsfont} +\usepackage[amsmath,thmmarks,thref,framed]{ntheorem} +\usepackage[dvips]{graphics} +\usepackage{epsfig} +\usepackage{calc} +\usepackage{tabls} +\usepackage{slashbox} +\usepackage{times} +\usepackage{multicol} +\usepackage{tabularx} +\usepackage{textcomp} + +\usepackage{pst-all} + +\usepackage[a4paper]{geometry} + + +\date{} +\geometry{hmargin=1cm, tmargin=1cm,bmargin=1.5cm} +\begin{document} +\title{UE MESI, Master IMR 2ème année.\\ + Novembre 2012 (durée 1h). J.-F. Couchot,} + +\maketitle +\vspace{-5em} +\begin{tabular}{ll} +Nom:& ........................................\\ +Prénom:& ........................................\\ +\end{tabular} + + +On s'intéresse à résoudre une équation de la forme $f(x)=0$ par la +méthode de Müller. Dans cette méthode, on considère que l'on +a le triplet de points $(x_{n-2},x_{n-1},x_{n})$. Pour calculer +$x_{n+1}$, on fait comme suit: +\begin{enumerate} +\item \label{itm:1} on approche $f(x)$ par un polynôme $P(x)$ aux points + $(x_{n-2},x_{n-1},x_{n})$, +\item\label{itm:2} on résout l'équation $P(x)=0$. La racine la plus proche de $x_n$ est $x_{n+1}$; +\item on recommence avec le triplet $(x_{n-1},x_{n},x_{n+1})$\ldots +\end{enumerate} + + + + + + +Vos réponses seront données directement ci-dessous. + +\begin{enumerate} + +\item En utilisant une base de Lagrange, montrer que le polynôme $P(x)$ + obtenu à l'étape 1. de la première itération est défini par +$$ +P(x) = \dfrac{(x-x_{n-1})(x-x_{n-2})f(x_n)}{(x_n-x_{n-1})(x_n-x_{n-2})} + +\dfrac{(x-x_{n})(x-x_{n-2})f(x_{n-1})}{(x_{n-1}-x_{n})(x_{n-1}-x_{n-2})} + +\dfrac{(x-x_{n})(x-x_{n-1})f(x_{n-2})}{(x_{n-2}-x_{n})(x_{n-2}-x_{n-1})} +$$ + +\vspace{4cm} + +\item Montrer que le polynôme de la question précédente est de degré 2. + Est-ce cohérent avec le fait qu'on veuille approximer $f$ en trois points? +\vspace{4cm} + +\item Montrer que le polynôme de la première question peut s'écrire + sous la forme $P(x) = a_n x^2 + b_n x + c $ où +\begin{eqnarray*} +a_n & = & \dfrac{f(x_n)}{(x_n-x_{n-1})(x_n-x_{n-2})} + +\dfrac{f(x_{n-1})}{(x_{n-1}-x_{n})(x_{n-1}-x_{n-2})} + +\dfrac{f(x_{n-2})}{(x_{n-2}-x_{n})(x_{n-2}-x_{n-1})} \\ +b_n & = &-\dfrac{f(x_n)(x_{n-1}+x_{n-2})}{(x_n-x_{n-1})(x_n-x_{n-2})} - +\dfrac{f(x_{n-1})(x_{n}+x_{n-2})}{(x_{n-1}-x_{n})(x_{n-1}-x_{n-2})} - +\dfrac{f(x_{n-2})(x_{n}+x_{n-1})}{(x_{n-2}-x_{n})(x_{n-2}-x_{n-1})} \\ +c_n & = & \dfrac{f(x_n)x_{n-1}x_{n-2}}{(x_n-x_{n-1})(x_n-x_{n-2})} + +\dfrac{f(x_{n-1})x_{n}x_{n-2}}{(x_{n-1}-x_{n})(x_{n-1}-x_{n-2})} + +\dfrac{f(x_{n-2})x_{n}x_{n-1}}{(x_{n-2}-x_{n})(x_{n-2}-x_{n-1})} +\end{eqnarray*} +\vspace{8cm} + + +\item Exprimer les deux racines $x'_{n}$ et $x''_{n}$ du polynôme précédent +en fonctions de $a_n$, $b_n$ et $c_n$ lorsqu'on itère dans les réels. +\vspace{3cm} + +\item Comment est alors défini $x_{n+1}$? +\vspace{3cm} + +\item On pourrait montrer que l'ordre de la convergence est 1,84. Comparer cette vitesse de convergence avec celle de Newton et celle de Lagrange. +\vspace{3cm} + + +\item Donner le code Python de la fonction + $\verb+[n,X] = iteration_muller(x+_{\verb+0+},\verb+x+_{\verb+1+},\verb+x+_{\verb+2+}\verb+,m,epsilon,f)+$ où +\begin{itemize} +\item $\verb+x+_{\verb+0+}$, $\verb+x+_{\verb+1+}$ et $\verb+x+_{\verb+2+}$ sont les trois premières valeurs des itérés, \verb+m+ + est le nombre maximal + d'itérations, \texttt{epsilon} est la précision souhaitée + et \verb+f+ la fonction à itérer; +\item \verb+n+ est le nombre d'itérations réalisées pour que +\verb+f(+$\verb+x+_{\verb+n+}$\verb+)+=0 ou que +$|\verb+x+_{\verb+n+}- \verb+x+_{\verb+n-1+}| \leq \verb+epsilon+$, \verb+n+ étant inférieur à \verb+m+ et \verb+X+ est + le vecteur contenant les + valeurs $\verb+x+_{\verb+0+},\ldots,\verb+x+_{\verb+n+}$. +\end{itemize} +\end{enumerate} + + +\end{document} \ No newline at end of file diff --git "a/tel/IntroAuxM\303\251thodesNum\303\251riques.pdf" "b/tel/IntroAuxM\303\251thodesNum\303\251riques.pdf" new file mode 100644 index 0000000..26f88e8 Binary files /dev/null and "b/tel/IntroAuxM\303\251thodesNum\303\251riques.pdf" differ diff --git a/tel/TPmatlab/equation_differentielle/TP5a/demo_affiche_courbe.m b/tel/TPmatlab/equation_differentielle/TP5a/demo_affiche_courbe.m new file mode 100644 index 0000000..ff73a34 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5a/demo_affiche_courbe.m @@ -0,0 +1,53 @@ +% demo_affiche_courbe : ce script permet d'afficher des courbes +% correspondant à des approximations d'équations différentielles. +% +% paramètres d'entrée : +% * nombre de courbes ; +% * la borne supérieure de l'intervalle d'étude T ; +% * pour chacune de ces courbes, X un vecteur contenant les valeurs, déja calculées, +% aux instants 0,h,2h,...,Nh=T. +% + +%nombre de points pour l'affichage défini par ngra +ngra=1000; + +T=input('entrez la borne supérieur de l''intervalle : '); +nc=input('entrez le nombre de courbes d''erreurs à tracer : '); +YY=cell(nc,1); +xx=cell(nc,1); +yy=cell(nc,1); +for k=1:nc + YY{k}=input(['entrez le vecteur des valeurs de la courbe numero ',num2str(k),': ']); +end +for k=1:nc + N=length(YY{k})-1; + if (N=nmax)|(abs(v-oldv)<=epsilon)); + end + yapp=v; + if (n>=nmax) + error('attention, n>nmax pour la méthode de Newton : arrêt du programme'); + end + Y(i+1)=yapp; +end + + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5a/resout_runge_kutta2.m b/tel/TPmatlab/equation_differentielle/TP5a/resout_runge_kutta2.m new file mode 100644 index 0000000..c2d02c9 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5a/resout_runge_kutta2.m @@ -0,0 +1,57 @@ +function Y=resout_runge_kutta2(N,T,y0,fcn) +% +% +% +% resout_runge_kutta2 : résolution d'une équation différentielle par la méthode de Runge Kutta 2 (explicite). +% +% ********************************************************* +% +% Y=resout_runge_kutta2(N,T,y0,fcn) calcule les valeurs aux instants 0,h,2h,...,Nh=T par la méthode +% de Runge Kutta d'ordre 2 (explicite) pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0. +% +% * variables d'entrées : +% N : nombre de valeurs calculées ; +% T : borne supérieure de l'intervalle ; +% y0 : valeur de y à t=0 ; +% fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f. +% +% * variables de sortie : +% Y : valeurs calculée aux instants 0,h,2h,...,Nh=T. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + +% Contrôles d'entrée + +% nombre d'arguments +if nargin~=4 + error('nombre d''arguments de la fonction incorrect'); +end + +% Corps d'algorithme +h=T/N; +Y=zeros(1,N+1); +Y(1)=y0; +auxiY=y0; +for i=1:N + t=h*(i-1); + y=auxiY; + k1=h*eval(fcn); + t=t+h; + y=auxiY+k1; + k2=h*eval(fcn); + auxiY=auxiY+(k1+k2)/2; + Y(i+1)=auxiY; +end + + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5a/resout_runge_kutta4.m b/tel/TPmatlab/equation_differentielle/TP5a/resout_runge_kutta4.m new file mode 100644 index 0000000..e4b2874 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5a/resout_runge_kutta4.m @@ -0,0 +1,62 @@ +function Y=resout_runge_kutta4(N,T,y0,fcn) +% +% +% +% resout_runge_kutta4 : résolution d'une équation différentielle par la méthode de Runge Kutta 4 (explicite). +% +% ********************************************************* +% +% Y=resout_runge_kutta4(N,T,y0,fcn) calcule les valeurs aux instants 0,h,2h,...,Nh=T par la méthode +% de Runge Kutta d'ordre 4 (explicite) pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0. +% +% * variables d'entrées : +% N : nombre de valeurs calculées ; +% T : borne supérieure de l'intervalle ; +% y0 : valeur de y à t=0 ; +% fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f. +% +% * variables de sortie : +% Y : valeurs calculée aux instants 0,h,2h,...,Nh=T. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + +% Contrôles d'entrée + +% nombre d'arguments +if nargin~=4 + error('nombre d''arguments de la fonction incorrect'); +end + +% Corps d'algorithme +h=T/N; +Y=zeros(1,N+1); +Y(1)=y0; +auxiY=y0; +for i=1:N + t=h*(i-1); + y=auxiY; + k1=h*eval(fcn); + t=t+h/2; + y=auxiY+k1/2; + k2=h*eval(fcn); + y=auxiY+k2/2; + k3=h*eval(fcn); + t=h*i; + y=auxiY+k3; + k4=h*eval(fcn); + auxiY=auxiY+(k1+2*k2+2*k3+k4)/6; + Y(i+1)=auxiY; +end + + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5a/resout_theta_methode.m b/tel/TPmatlab/equation_differentielle/TP5a/resout_theta_methode.m new file mode 100644 index 0000000..36de271 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5a/resout_theta_methode.m @@ -0,0 +1,81 @@ +function Y=resout_theta_methode(theta,N,T,y0,fcn,nmax,epsilon) +% +% +% +% resout_theta_methode : résolution d'une équation différentielle par la theta-méthode. +% +% ********************************************************* +% +% Y=resout_theta_methode(theta,N,T,y0,fcn,nmax,epsilon) calcule les valeurs aux instants 0,h,2h,...,Nh=T par +% la theta-méthode +% pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0. +% A chaque itération, l'équation y_{n+1}=y_n+theta hf(t_{n+1},y_{n+1})+(1-theta)hf(t_{n},y_{n}) +% est résolue par la méthode de Newton +% où la valeur initiale est choisie égale à y_n. +% +% La dérivée partielle df/dy est calculée en symbolique. +% +% * variables d'entrées : +% theta est un réel dans [0,1] ; +% N : nombre de valeurs calculées ; +% T : borne supérieure de l'intervalle ; +% y0 : valeur de y à t=0 ; +% fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f +% epsilon : précision souhaitée (pour la méthode de Newton) +% nmax : nombre maximal d'itérations (pour la méthode de Newton) +% +% * variables de sortie : +% Y : valeurs calculée aux instants 0,h,2h,...,Nh=T. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +% Contrôles d'entrée + +% nombre d'arguments +if nargin~=7 + error('nombre d''arguments de la fonction incorrect'); +end + +% Corps d'algorithme +h=T/N; +Y=zeros(1,N+1); +Y(1)=y0; +auxiY=y0; +syms t y; +g=eval(fcn); +fcnpsymb=diff(g,y); +fcnp=char(fcnpsymb); +yapp=y0; +for i=1:N + y=yapp; + t=h*(i-1); + auxidud=h*(1-theta)*eval(fcn); + t=h*i; + fx=-h*theta*eval(fcn)-auxidud; + fxp=1-h*theta*eval(fcnp); + test=((fx)~=0); + n=0; + v=yapp; + while(test) + oldv=v; + v=v-fx/fxp; + y=v; + t=h*i; + fx=v-yapp-h*theta*eval(fcn)-auxidud; + fxp=1-h*theta*eval(fcnp); + n=n+1; + test=~((fx==0)|(n>=nmax)|(abs(v-oldv)<=epsilon)); + end + yapp=v; + if (n>=nmax) + error('attention, n>nmax pour la méthode de Newton : arrêt du programme'); + end + Y(i+1)=yapp; +end diff --git a/tel/TPmatlab/equation_differentielle/TP5b/demo_log_erreur_schema.m b/tel/TPmatlab/equation_differentielle/TP5b/demo_log_erreur_schema.m new file mode 100644 index 0000000..f919e69 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5b/demo_log_erreur_schema.m @@ -0,0 +1,100 @@ +% demo_log_erreur_schema : ce script permet de déterminer +% l'ordre du schéma de façon graphique +% et numérique des cinq méthodes suivantes : +% * Euler explicite ; +% * Euler implicite ; +% * Runge Kutta explicite d'ordre 2 +% * theta_méthode +% * Runge Kutta explicite d'ordre 4 +% +% paramètres d'entrée : +% nmin et nmax : deux entiers entre lesquels varie n ; +% T : borne supérieure de l'intervalle ; +% y0 : valeur de y à t=0 ; +% theta est un réel dans [0,1] ; +% epsilon : précision souhaitée (pour la méthode de Newton, pour les méthodes implicites) +% nmax : nombre maximal d'itérations (pour la méthode de Newton, pour les méthodes implicites) +% fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f + +pmin=input('entrez l''exposant pmin tel que le plus petit entier de calcul soit 10^pmin : '); +pmax=input('entrez l''exposant pmax tel que le plus grand entier de calcul soit 10^pmax : '); +nc=input('entrez le nombre de calcul : '); +T=input('entrez la borne supérieur de l''intervalle : '); +y0=input('entrez la condition initiale : '); +theta=input('entrez le paramètre theta : '); +epsilon=input('entrez epsilon (pour la méthode de Newton, pour les méthodes implicites) : '); +nmax=input('entrez le nombre maximum d''itérations (pour la méthode de Newton, pour les méthodes implicites) : '); +fcn=input('entrez fcn (sous la forme ''f(t,y)'') : '); + +syms t y; +g=eval(fcn); +fcnpsymb=diff(g,y); +fcnp=char(fcnpsymb); +x=[]; +y1=[]; +y2=[]; +y3=[]; +y4=[]; +y5=[]; +for N=logspace(pmin,pmax,nc) + NN=fix(N); + x=[x log10(NN)]; + Ya1=resout_euler_explicite(NN,T,y0,fcn); + Yb1=resout_euler_explicite(2*NN,T,y0,fcn); + Ya2=resout_euler_implicite_bis(NN,T,y0,fcn,fcnp,nmax,epsilon); + Yb2=resout_euler_implicite_bis(2*NN,T,y0,fcn,fcnp,nmax,epsilon); + Ya3=resout_runge_kutta2(NN,T,y0,fcn); + Yb3=resout_runge_kutta2(2*NN,T,y0,fcn); + Ya4=resout_theta_methode_bis(theta,NN,T,y0,fcn,fcnp,nmax,epsilon); + Yb4=resout_theta_methode_bis(theta,2*NN,T,y0,fcn,fcnp,nmax,epsilon); + Ya5=resout_runge_kutta4(NN,T,y0,fcn); + Yb5=resout_runge_kutta4(2*NN,T,y0,fcn); + auxi1=-log10(max(abs(Ya1-Yb1(1:2:end)))); + auxi2=-log10(max(abs(Ya2-Yb2(1:2:end)))); + auxi3=-log10(max(abs(Ya3-Yb3(1:2:end)))); + auxi4=-log10(max(abs(Ya4-Yb4(1:2:end)))); + auxi5=-log10(max(abs(Ya5-Yb5(1:2:end)))); + y1=[y1 auxi1]; + y2=[y2 auxi2]; + y3=[y3 auxi3]; + y4=[y4 auxi4]; + y5=[y5 auxi5]; + disp(' '); + disp(['calcul pour N=',int2str(NN),' terminé']); +end +[a1,b,r1]=regression_lineaire(x,y1); +[a2,b,r2]=regression_lineaire(x,y2); +[a3,b,r3]=regression_lineaire(x,y3); +[a4,b,r4]=regression_lineaire(x,y4); +[a5,b,r5]=regression_lineaire(x,y5); +clf; +hold on; +plot(x,y1,x,y2,x,y3,x,y4,x,y5); +legend('Euler explicite','Euler implicite','runge Kutta 2','theta-méthode','runge Kutta 4',0); +title('détermination des ordres des schémas'); +hold off; +disp(' '); +disp('pour avoir une estimation des ordres, tapez sur une touche'); +pause; +disp(' '); +disp('Euler explicite'); +disp(['ordre=',num2str(a1)]); +disp(['correlation=',num2str(r1)]); +disp(' '); +disp('Euler implicite'); +disp(['ordre=',num2str(a2)]); +disp(['correlation=',num2str(r2)]); +disp(' '); +disp('runge Kutta 2'); +disp(['ordre=',num2str(a3)]); +disp(['correlation=',num2str(r3)]); +disp(' '); +disp('theta-méthode'); +disp(['ordre=',num2str(a4)]); +disp(['correlation=',num2str(r4)]); +disp(' '); +disp('runge Kutta 4'); +disp(['ordre=',num2str(a5)]); +disp(['correlation=',num2str(r5)]); + + diff --git a/tel/TPmatlab/equation_differentielle/TP5b/resout_euler_implicite_bis.m b/tel/TPmatlab/equation_differentielle/TP5b/resout_euler_implicite_bis.m new file mode 100644 index 0000000..058f68b --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5b/resout_euler_implicite_bis.m @@ -0,0 +1,76 @@ +function Y=resout_euler_implicite_bis(N,T,y0,fcn,fcnp,nmax,epsilon) +% +% +% +% resoud_euler_implicite_bis : résolution d'une équation différentielle par la méthode d'Euler implicite. +% +% ********************************************************* +% +% X=resout_euler_implicite_bis(N,T,y0,fcn,fcnp,nmax,epsilon) calcule les valeurs aux instants 0,h,2h,...,Nh=T +% par la méthode d'Euler implicite pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0. +% A chaque itération, l'équation y_{n+1}=y_n+hf(t_{n+1},y_{n+1}) est résolue par la méthode de Newton +% où la valeur initiale est choisie égale à y_n. +% +% La dérivée partielle df/dy est entrée comme paramètre. +% +% * variables d'entrées : +% N : nombre de valeurs calculées ; +% T : borne supérieure de l'intervalle ; +% y0 : valeur de y à t=0 ; +% fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction +% fcnp : chaine de caractères (de type 'f(t,y)') représentant la fonction df/dy +% epsilon : précision souhaitée (pour la méthode de Newton) +% nmax : nombre maximal d'itérations (pour la méthode de Newton) +% +% * variables de sortie : +% Y : valeurs calculée aux instants 0,h,2h,...,Nh=T. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +% Contrôles d'entrée + +% nombre d'arguments +if nargin~=7 + error('nombre d''arguments de la fonction incorrect'); +end + +% Corps d'algorithme +h=T/N; +Y=zeros(1,N+1); +Y(1)=y0; +auxiY=y0; +yapp=y0; +for i=1:N + t=h*i; + y=yapp; + fx=-h*eval(fcn); + fxp=1-h*eval(fcnp); + test=((fx)~=0); + n=0; + v=yapp; + while(test) + oldv=v; + v=v-fx/fxp; + y=v; + fx=v-yapp-h*eval(fcn); + fxp=1-h*eval(fcnp); + n=n+1; + test=~((fx==0)|(n>=nmax)|(abs(v-oldv)<=epsilon)); + end + yapp=v; + if (n>=nmax) + error('attention, n>nmax pour la méthode de Newton : arrêt du programme'); + end + Y(i+1)=yapp; +end + + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5b/resout_theta_methode_bis.m b/tel/TPmatlab/equation_differentielle/TP5b/resout_theta_methode_bis.m new file mode 100644 index 0000000..c65ee7c --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5b/resout_theta_methode_bis.m @@ -0,0 +1,78 @@ +function Y=resout_theta_methode_bis(theta,N,T,y0,fcn,fcnp,nmax,epsilon) +% +% +% +% resout_theta_methode_bis : résolution d'une équation différentielle par la theta-méthode. +% +% ********************************************************* +% +% Y=resout_theta_methode_bis(theta,N,T,y0,fcn,nmax,epsilon) calcule les valeurs aux instants 0,h,2h,...,Nh=T par +% la theta-méthode +% pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0. +% A chaque itération, l'équation y_{n+1}=y_n+theta hf(t_{n+1},y_{n+1})+(1-theta)hf(t_{n},y_{n}) +% est résolue par la méthode de Newton +% où la valeur initiale est choisie égale à y_n. +% +% La dérivée partielle df/dy est entrée comme paramètre. +% +% * variables d'entrées : +% theta est un réel dans [0,1] ; +% N : nombre de valeurs calculées ; +% T : borne supérieure de l'intervalle ; +% y0 : valeur de y à t=0 ; +% fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f +% fcnp : chaine de caractères (de type 'f(t,y)') représentant la fonction df/dy +% epsilon : précision souhaitée (pour la méthode de Newton) +% nmax : nombre maximal d'itérations (pour la méthode de Newton) +% +% * variables de sortie : +% Y : valeurs calculée aux instants 0,h,2h,...,Nh=T. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +% Contrôles d'entrée + +% nombre d'arguments +if nargin~=8 + error('nombre d''arguments de la fonction incorrect'); +end + +% Corps d'algorithme +h=T/N; +Y=zeros(1,N+1); +Y(1)=y0; +auxiY=y0; +yapp=y0; +for i=1:N + y=yapp; + t=h*(i-1); + auxidud=h*(1-theta)*eval(fcn); + t=h*i; + fx=-h*theta*eval(fcn)-auxidud; + fxp=1-h*theta*eval(fcnp); + test=((fx)~=0); + n=0; + v=yapp; + while(test) + oldv=v; + v=v-fx/fxp; + y=v; + t=h*i; + fx=v-yapp-h*theta*eval(fcn)-auxidud; + fxp=1-h*theta*eval(fcnp); + n=n+1; + test=~((fx==0)|(n>=nmax)|(abs(v-oldv)<=epsilon)); + end + yapp=v; + if (n>=nmax) + error('attention, n>nmax pour la méthode de Newton : arrêt du programme'); + end + Y(i+1)=yapp; +end diff --git a/tel/TPmatlab/equation_differentielle/TP5c/calcul_developpement_limite.m b/tel/TPmatlab/equation_differentielle/TP5c/calcul_developpement_limite.m new file mode 100644 index 0000000..05a50d5 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5c/calcul_developpement_limite.m @@ -0,0 +1,38 @@ +function dly=calcul_developpement_limite(t,y,h,p); + +% +% calcul_developpement_limite : calcul symbolique du développement limité ... +% +% ********************************************************* +% +% dly=calcul_developpement_limite(t,y,h,p) est le développement limité +% formel de y, solution de y'(t)=f(t,y(t)) autour de y() à l'ordre p. +% variables d'entrées : +% * entier naturel p +% * symbolique t , y et h +% variables de sortie : dly=alcul_developpement_limite(t,y,h,p) +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +syms yb tb hb; +auxifacto=sym(1); +factoriel=[]; +for j=1:p + auxifacto=auxifacto*j; + factoriel=[factoriel 1/auxifacto]; +end +puissanceh=hb.^(1:p); +F=calcul_fonction_fp(tb,yb,p-1); +if (p==0) + dly=yb; +else + dly=yb+sum(factoriel.*puissanceh.*F); +end +dly=subs(dly,{'tb','yb','hb'},{char(t),char(y),char(h)}); diff --git a/tel/TPmatlab/equation_differentielle/TP5c/calcul_fonction_fp.m b/tel/TPmatlab/equation_differentielle/TP5c/calcul_fonction_fp.m new file mode 100644 index 0000000..85d4c52 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5c/calcul_fonction_fp.m @@ -0,0 +1,50 @@ +function F=calcul_fonction_fp(t,y,p); + +% +% calcul_fonction_fp : définition symbolique des fonctions f_0,...,f_p +% +% ********************************************************* +% +% F=calcul_fonction_fp(t,y,p) +% calcule les fonctions f_0, f_1, ..., f_p définies par récurrence par +% f_0(t,y)=f(t,y) +% pour k dans {0,...,p-1}, +% f_{k+1}=diff(f_k,t)(t,y)+diff(f_k,y)(t,y)*f(t,y). +% Chaque foncton f_k est stockée dans F(k+1) pour k dans {0,...,p-1}. +% * variables d'entrée : +% * entier naturel p +% * symbolique t , y et h +% * variables de sortie : F=[f_0 f_1 ... f_p] +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +if (p>=0) + syms yb tb; + fcn=sym('f(tb,yb)'); + fauxi=fcn; + F=zeros(1,p+1); + F=sym(F); + F(1)=fcn; + for k=0:p-1 + fauxi=simplify(diff(fauxi,tb)+diff(fauxi,yb)*fcn); + F(k+2)=fauxi; + end + for k=1:p+1 + F(k)=subs(F(k),{'tb','yb'},{char(t),char(y)});; + end +else + F=sym([]); +end + + + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5c/demo_determine_ordre.m b/tel/TPmatlab/equation_differentielle/TP5c/demo_determine_ordre.m new file mode 100644 index 0000000..5929c32 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5c/demo_determine_ordre.m @@ -0,0 +1,34 @@ +% ce script determine_ordre calcule les erreurs de consistance pour les +% schémas numériques d'Euler explicite, Runge Kutta 2 et 4. + +disp('calcul en cours ....'); +syms t y h; +dly=calcul_developpement_limite(t,y,h,6); +phiee=phi_euler_explicite_symbolique(t,y,h); +phirk2=phi_runge_kutta2_symbolique(t,y,h); +phirk4=phi_runge_kutta4_symbolique(t,y,h); +ynpuee=y+h*phiee; +ynpurk2=y+h*phirk2; +ynpurk4=y+h*phirk4; +see=simplify(taylor(dly,h,3)-taylor(ynpuee,h,3)); +srk2=simplify(taylor(dly,h,4)-taylor(ynpurk2,h,4)); +srk4=simplify(taylor(dly,h,6)-taylor(ynpurk4,h,6)); +disp('développement limité à l''ordre 5 : '); +pretty(expand(taylor(dly,h,6))); +disp('appuyez sur une touche pour continuer.'); +pause; +disp('erreur de consistance pour Euler explicite (à l''ordre 2) : '); +pretty(see); +disp('appuyez sur une touche pour continuer.'); +pause; +disp('erreur de consistance pour RK2 (à l''ordre 3) : '); +pretty(srk2); +disp('appuyez sur une touche pour continuer.'); +pause; +disp('erreur de consistance pour RK4 (à l''ordre 5) : '); +pretty(srk4); +disp('appuyez sur une touche pour continuer.'); +pause; + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5c/phi_euler_explicite_symbolique.m b/tel/TPmatlab/equation_differentielle/TP5c/phi_euler_explicite_symbolique.m new file mode 100644 index 0000000..940b779 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5c/phi_euler_explicite_symbolique.m @@ -0,0 +1,29 @@ +function phi=phi_euler_explicite_symbolique(t,y,h); + +% +% phi_euler_explicite_symbolique : définition symbolique du schéma d'Euler (Phi) +% +% ********************************************************* +% +% phi=phi_euler_explicite_symbolique(t,y,h) +% calcule la fonction (en symbolique) : (t,y)-> Phi(t,y) telle +% le schéma d'Euler explicite s'écrive sous la forme : +% y_{i+1}=y_{i}+h Phi(t_i,y_i,h) +% variables d'entrée : t , y et h (symbolique)) +%% variables de sortie : res=Phi(t,y,h) +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +syms yb tb hb; +fcn=sym('f(tb,yb)'); +phi=subs(fcn,{'tb','yb'},{char(t),char(y)}); + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5c/phi_runge_kutta2_symbolique.m b/tel/TPmatlab/equation_differentielle/TP5c/phi_runge_kutta2_symbolique.m new file mode 100644 index 0000000..e81b984 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5c/phi_runge_kutta2_symbolique.m @@ -0,0 +1,37 @@ +function phi=phi_runge_kutta2_symbolique(t,y,h); + +% +% phi_runge_kutta2_symbolique : définition symbolique de Runge Kutta 2 (Phi) +% +% ********************************************************* +% +% phi=phi_runge_kutta2_symbolique(t,y,h) +% calcule la fonction (en symbolique) : (t,y)-> Phi(t,y) telle +% le schéma de Runge Kutta 2 s'écrive sous la forme : +% y_{i+1}=y_{i}+h Phi(t_i,y_i,h) +% variables d'entrée : t , y et h (symbolique)) +%% variables de sortie : res=Phi(t,y,h) +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +syms yb tb hb; +fcn=sym('f(tb,yb)'); +k1tilde=fcn; +fcnb=subs(fcn,{'tb','yb'},{tb+hb,yb+hb*k1tilde}); +k2tilde=fcnb; +phi=(k1tilde+k2tilde)/sym(2); +phi=subs(phi,{'tb','yb','hb'},{char(t),char(y),char(h)}); + + + + + + + diff --git a/tel/TPmatlab/equation_differentielle/TP5c/phi_runge_kutta4_symbolique.m b/tel/TPmatlab/equation_differentielle/TP5c/phi_runge_kutta4_symbolique.m new file mode 100644 index 0000000..3bce3b4 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5c/phi_runge_kutta4_symbolique.m @@ -0,0 +1,34 @@ +function phi=phi_runge_kutta4_symbolique(t,y,h); + +% +% phi_runge_kutta4_symbolique : définition symbolique de Runge Kutta 4 (Phi) +% +% ********************************************************* +% +% phi=phi_runge_kutta4_symbolique(t,y,h) +% calcule la fonction (en symbolique) : (t,y)-> Phi(t,y) telle +% le schéma de Runge Kutta 4 s'écrive sous la forme : +% y_{i+1}=y_{i}+h Phi(t_i,y_i,h) +% variables d'entrée : t , y et h (symbolique)) +%% variables de sortie : res=Phi(t,y,h) +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +syms yb tb hb; +fcn=sym('f(tb,yb)'); +k1tilde=fcn; +fcnb=subs(fcn,{'tb','yb'},{tb+1/2*hb,yb+1/2*hb*k1tilde}); +k2tilde=fcnb; +fcnb=subs(fcn,{'tb','yb'},{tb+1/2*hb,yb+1/2*hb*k2tilde}); +k3tilde=fcnb; +fcnb=subs(fcn,{'tb','yb'},{tb+hb,yb+hb*k3tilde}); +k4tilde=fcnb; +phi=1/6*(k1tilde+2*k2tilde+2*k3tilde+k4tilde); +phi=subs(phi,{'tb','yb','hb'},{char(t),char(y),char(h)}); diff --git a/tel/TPmatlab/equation_differentielle/TP5d/demo_trace_systeme_mecanique.m b/tel/TPmatlab/equation_differentielle/TP5d/demo_trace_systeme_mecanique.m new file mode 100644 index 0000000..5975791 --- /dev/null +++ b/tel/TPmatlab/equation_differentielle/TP5d/demo_trace_systeme_mecanique.m @@ -0,0 +1,106 @@ +% ce script demo_trace_systeme_mecanique calcule les +% solutions du système à deux degrés de liberté régi par le système : +% x1''(t)+(k1+k2)x1-k2x2+c1x'1+d1(x'1)^3=f1 +% x2''(t)-k2x1+(k2+k3)x2+c2x'2+d2(x'2)^3=f2 +% déterminée par les schémas d'Euler implicite et de Runge Kutta 4. +% +% +% * les paramètres d'entrée sont : +% k : le vecteur contenant k1, k2 et k3, +% c : le vecteur contenant c1 et c2, +% d : le vecteur contenant d1, d2, +% y0 : le vecteur contenant les conditions initiales (x1(0), x1'(0), x2(0), x2'(0)) +% fcn1 : chaine de caractères (de type 'f1(t)') représentant la fonction f1, +% fcn2 : chaine de caractères (de type 'f2(t)') représentant la fonction f2, +% N : un entier non nul, +% T : la borne supérieur de l'intervalle d'étude. +% +% * le script affiche les approximation de x1, x1', x2 et x2' sur [0,T]. + +clear all; +%nombre de points pour l'affichage défini par ngra +ngra=1000; + +% entree des paramètres +k=input('entrez le vecteur [k1,k2,k3] : '); +c=input('entrez le vecteur [c1,c2]: '); +d=input('entrez le vecteur [d1,d2]: '); +y0=input('entrez le vecteur des conditions initiales [x1(0),x1''(0),x2(0),x2''(0)]: '); +fcn1=input('entrez la fonction f1 sous la forme ''f1(t)'': '); +fcn2=input('entrez la fonction f2 sous la forme ''f2(t)'': '); +N=input('entrez le nombre de points de discretisation : '); +T=input('entrez T : '); + + +% définition de la fonction F telle que le système soit équivalent à Y'(t)=F(t,Y(t)) dans R^4. +% F est defini en fonction de t et de y, vecteur à 4 composantes. +A=-[0 -1 0 0; + k(1)+k(2) c(1) -k(2) 0; + 0 0 0 -1; + -k(2) 0 k(2)+k(3) c(2)]; +H='-[0;d(1)*(y(2))^3;0;d(2)*(y(4))^3]'; +G=['[0; ',fcn1,';','0;',fcn2,']']; +fcn=[H,'+',G,'+A*y']; +h=T/N; + +% calcul de la solution par Euler explicite. +Yeuler=zeros(4,N+1); +Yeuler(:,1)=y0'; +auxiY=y0'; +for i=1:N + t=h*(i-1); + y=auxiY; + auxiY=auxiY+h*eval(fcn); + Yeuler(:,i+1)=auxiY; +end + +% calcul de la solution par Runge Kutta. +Yrungekutta=zeros(4,N+1); +Yrungekutta(:,1)=y0'; +auxiY=y0'; +for i=1:N + t=h*(i-1); + y=auxiY; + k1=h*eval(fcn); + t=t+h/2; + y=auxiY+k1/2; + k2=h*eval(fcn); + y=auxiY+k2/2; + k3=h*eval(fcn); + t=h*i; + y=auxiY+k3; + k4=h*eval(fcn); + auxiY=auxiY+(k1+2*k2+2*k3+k4)/6; + Yrungekutta(:,i+1)=auxiY; +end + + +% préparation des vecteurs de tracé. +if (N=nmax +% * X est une vecteur qui contient les valeurs de x_0 à x_n +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + + + + +% Contrôles d'entrée +if (a>=b) + error('il faut a=nmax)|(abs(x-xm1)<=epsilon)); +end +if (n>=nmax) + disp('attention, n>nmax'); +end + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_dichotomie.m b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_dichotomie.m new file mode 100644 index 0000000..1a009cd --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_dichotomie.m @@ -0,0 +1,60 @@ +function [n,X]=iteration_dichotomie(a,b,nmax,epsilon,fcn) + + +% iteration_dichotomie : calcul d'une racine d'une fonction par dichotomie +% +% ********************************************************* +% +% [n,X]=iteration_dichotomie(a,b,nmax,epsilon,fcn) renvoie +% les itérés de la méthode de dichotomie +% +% variables d'entrées : +% * a,b : borne de l'intervalle (avec f(a)f(b)<=0) +% * nmax : nombre maximal d'itérations +% * epsilon : précision souhaitée +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% +% variables de sortie : +% * n est l'indice tel que (1/2^n)(b-a)=b) + error('il faut a0) + error('il faut f(a)f(b)<=0'); +end + +% corps d'algorithme +X=(a+b)/2; +n=fix(-log(epsilon/(b-a))/log(2))+1; +if (n>nmax) + error('il faut n plus petit que nmax'); +end +for k=1:n + x=(a+b)/2; + X=[X x]; + auxif=feval(fcn,a)*feval(fcn,x); + if (auxif<=0) + b=x; + else + a=x; + end +end diff --git a/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_newton.m b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_newton.m new file mode 100644 index 0000000..f53041a --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_newton.m @@ -0,0 +1,71 @@ +function [n,X]=iteration_newton(x0,nmax,epsilon,fcn) + + +% iteration_newton : calcul d'une racine la méthode de Newton +% +% ********************************************************* +% +% [n,X]=iteration_newton(x0,nmax,epsilon,fcn) renvoie +% les itérés de la méthode de Newton +% (la dérivée f' est calculée symboliquement) +% +% variables d'entrées : +% * nmax : nombre maximal d'itérations +% * x0 : le premier terme +% * epsilon : précision souhaitée +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% +% variables de sortie : +% * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|=nmax +% * X est une vecteur qui contient les valeurs de x_0 à x_n +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + +% Contrôles d'entrée +qs=exist(fcn); +if (qs~=2) & (qs~=5) & (qs~=1) + error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)'); +end + +% corps d'algorithme +syms vx +g=feval(fcn,vx); +fcnpsymb=diff(g,vx); +n=0; +test=(feval(fcn,x0)~=0); +x=x0; +fx=feval(fcn,x); +fxp=subs(fcnpsymb,vx,x); +if ~(isnumeric(fxp)); + fxp=eval(fxp); +end +X=x0; +while(test) + xnp1=x-fx/fxp; + n=n+1; + fxnp1=feval(fcn,xnp1); + test=~((fxnp1==0)|(n>=nmax)|(abs(x-xnp1)<=epsilon)); + x=xnp1; + fx=feval(fcn,x); + fxp=subs(fcnpsymb,vx,x); + if ~(isnumeric(fxp)); + fxp=eval(fxp); + end + X=[X x]; +end +if (n>=nmax) + disp('attention, n>nmax'); +end + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_regula_falsi.m b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_regula_falsi.m new file mode 100644 index 0000000..d299fee --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_regula_falsi.m @@ -0,0 +1,74 @@ +function [n,X]=iteration_regula_falsi(xm1,x0,nmax,epsilon,fcn) + + +% iteration_regula_falsi : calcul d'une racine par regula falsi +% +% ********************************************************* +% +% [n,X]=iteration_regula_falsi(x0,xm1,nmax,epsilon,fcn) renvoie +% les itérés de la méthode de la la fausse position +% +% variables d'entrées : +% * nmax : nombre maximal d'itérations +% * xm1 et x0 : les deux premiers termes (x_{-1} et x_0) +% * epsilon : précision souhaitée +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% +% variables de sortie : +% * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|=nmax +% * X est une vecteur qui contient les valeurs de x_0 à x_n +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + +% Contrôles d'entrée +qs=exist(fcn); +if (feval(fcn,x0)*feval(fcn,xm1)>0) + error('il faut f(x0)f(x_(-1))<0'); +end +if (qs~=2) & (qs~=5) & (qs~=1) + error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)'); +end + +% corps d'algorithme +p=-1; +n=0; +test=(feval(fcn,x0)~=0); +x=x0; +xp=xm1; +fxp=feval(fcn,xm1); +fx=feval(fcn,x); +xmm1=xm1; +X=x0; +while(test) + xnp1=x-((x-xp)/(fx-fxp))*fx; + fxnp1=feval(fcn,xnp1); + test=~((fxnp1==0)|(n>=nmax)|(abs(x-xmm1)<=epsilon)); + xmm1=x; + x=xnp1; + fx=feval(fcn,x); + X=[X x]; + if (test) + if (feval(fcn,xmm1)*fx<0) + p=n; + xp=xmm1; + fxp=feval(fcn,xp); + end + end + n=n+1; +end +if (n>=nmax) + disp('attention, n>nmax'); +end + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_secante.m b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_secante.m new file mode 100644 index 0000000..fa67990 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4a/iteration_secante.m @@ -0,0 +1,64 @@ +function [n,X]=iteration_secante(xm1,x0,nmax,epsilon,fcn) + + +% iteration_secante : calcul d'une racine la méthode de la sécante +% +% ********************************************************* +% +% [n,X]=iteration_secante(x0,xm1,nmax,epsilon,fcn) renvoie +% les itérés de la méthode de la sécante +% +% variables d'entrées : +% * nmax : nombre maximal d'itérations +% * xm1 et x0 : les deux premiers termes (x_{-1} et x_0) +% * epsilon : précision souhaitée +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% +% variables de sortie : +% * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|=nmax +% * X est une vecteur qui contient les valeurs de x_0 à x_n +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + + +% Contrôles d'entrée +qs=exist(fcn); +if (qs~=2) & (qs~=5) & (qs~=1) + error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)'); +end + +% corps d'algorithme +n=0; +test=(feval(fcn,x0)~=0); +x=x0; +xmm1=xm1; +fxmm1=feval(fcn,xm1); +fx=feval(fcn,x); +X=x0; +while(test) + xnp1=x-((x-xmm1)/(fx-fxmm1))*fx; + n=n+1; + fxnp1=feval(fcn,xnp1); + test=~((fxnp1==0)|(n>=nmax)|(abs(x-xmm1)<=epsilon)); + xmm1=x; + fxmm1=fx; + x=xnp1; + fx=feval(fcn,x); + X=[X x]; +end +if (n>=nmax) + disp('attention, n>nmax'); +end + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4a/log_erreur_iteration.m b/tel/TPmatlab/equation_nonlineaire/TP4a/log_erreur_iteration.m new file mode 100644 index 0000000..7b6c603 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4a/log_erreur_iteration.m @@ -0,0 +1,39 @@ +function [ordre,corel]=log_erreur_iteration(l,X) + + +% log_erreur_iteration : détermination d'un ordre de convergence (équation non linéaire) +% +% ********************************************************* +% +% [ordre,corel]=log_erreur_iteration(l,X) renvoie +% l'ordre et la correlation à partir d'itérés déjà calculé par une méthode numérique de +% résolution de l'équation (scalaire) f(x)=0. +% +% variables d'entrées : +% * X : le vecteur des différentes valeurs de la suite x_n +% (déjà calculé par une méthode type dichotomie, Newton ....) +% * l : la limite présumée (souvent choisie égale à la dernière valeur de X) +% +% variables de sortie : +% * ordre est l'ordre mesuré (par droite de moindres carrés, à partir du nuage log-log) +% * corel : correlation mesurée +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% regression_lineaire +% +% ********************************************************* +% + +% corps d'algorithme +auxi=abs(X-l); +ind=find(auxi~=0); +auxi=auxi(ind); +auxi=log10(auxi); +x=auxi(1:end-1); +y=auxi(2:end); +[ordre,b,corel]=regression_lineaire(x,y); + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4b/demo_trace_iteration_point_fixe.m b/tel/TPmatlab/equation_nonlineaire/TP4b/demo_trace_iteration_point_fixe.m new file mode 100644 index 0000000..49ab608 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4b/demo_trace_iteration_point_fixe.m @@ -0,0 +1,74 @@ +% demo_trace_iteration_point_fixe : tracé des différents itérés d'une méthode de point fixe. +% +% à entrer : * la fonction g (telle que x_{n+1}=g(x_n)), +% * la valeur initiale (numérique ou symbolique) +% * p l'indice de la dernière valeur calculée. +% affiche le colimaçon : +% * avec les valeurs successives avec animation par movie +% * ou en statique + +clear g x0 X F k p; +disp('entrée de la fonction g'); +g=saisiefonction; +x0=input('entrez la première valeur : '); +p=input('entrez le nombre de valeurs : '); +ch=input('voulez vous le colimaçon en statique (entrez s) ou en animation (entrez a) : ','s'); +if isnumeric(x0) + X=zeros(1,p+1); + ch0=' (numérique)'; + vx0=num2str(x0); +else + X=sym(zeros(1,p+1)); + ch0=' (symbolique)'; + vx0=char(x0); +end +auxi=x0; +X(1)=x0; +for k=1:p + auxi=feval(g,auxi); + X(k+1)=auxi; +end +p=length(X)-1; +if ~isnumeric(x0) + X=double(vpa(X)); +end +xmin=min(X); +xmax=max(X); +if (xmin==xmax) + error('attention, suite constante : arrêt du programme'); +end +clf; +x2=xmax+0.1*(xmax-xmin); +x1=xmin-0.1*(xmax-xmin); +xgraph=x1:(x2-x1)/200:x2; +abscisse=zeros(1,2*p+1); +ordonnee=zeros(1,2*p+1); +abscisse(1:2:2*p+1)=X; +abscisse(2:2:2*p)=X(1:end-1); +ordonnee(1)=xmin; +ordonnee(2:2:2*p)=X(2:end); +ordonnee(3:2:2*p+1)=X(2:end); +hold on; +axis([x1 x2 x1 x2]); +plot(xgraph,feval(g,xgraph),'r'); +plot([x1 x2],[x1 x2],'g'); +if (strcmp(ch,'s')) + c=exist(g); + if (c==2) | (c==5) + titre=[g,'(x)']; + else + titre=formula(g); + end + + title([int2str(p),' itérés de la fonction ',titre, ' pour x0=',vx0,ch0]); + plot(abscisse,ordonnee,'b'); + legend('y=g(x)','y=x','itérés',0); +else + for k=0:p-1 + plot(abscisse(1:2*k+3),ordonnee(1:2*k+3)); + plot(X(k+1),X(k+2),'o b'); + F(k+1)=getframe; + end + movie(F,1,6); +end +hold off; diff --git a/tel/TPmatlab/equation_nonlineaire/TP4b/iteration_point_fixe.m b/tel/TPmatlab/equation_nonlineaire/TP4b/iteration_point_fixe.m new file mode 100644 index 0000000..77f5c33 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4b/iteration_point_fixe.m @@ -0,0 +1,57 @@ +function [n,X]=iteration_point_fixe(x0,nmax,epsilon,fcn) + + +% iteration_point_fixe : calcul des itérés de la méthode de point fixe +% +% ********************************************************* +% +% [n,X]=iteration_point_fixe(x0,nmax,epsilon,fcn) renvoie +% les itérés de la méthode de point fixe +% +% variables d'entrées : +% * nmax : nombre maximal d'itérations +% * x0 : terme initial +% * epsilon : précision souhaitée +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% +% variables de sortie : +% * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|=nmax +% * X est une vecteur qui contient les valeurs de x_0 à x_n +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% *********************************************************% + + + + + + +% Contrôles d'entrée +qs=exist(fcn); +if (qs~=2) & (qs~=5) & (qs~=1) + error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)'); +end + +% corps d'algorithme +n=0; +test=(feval(fcn,x0)-x0~=0); +x=x0; +X=x0; +xm1=x0; +while(test) + x=feval(fcn,x); + X=[X x]; + n=n+1; + test=~((feval(fcn,x)-x0==0)|(n>=nmax)|(abs(x-xm1)<=epsilon)); + xm1=x; +end +if (n>=nmax) + disp('attention, n>nmax'); +end + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4c/etude_fonction_particulier.m b/tel/TPmatlab/equation_nonlineaire/TP4c/etude_fonction_particulier.m new file mode 100644 index 0000000..72271f5 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4c/etude_fonction_particulier.m @@ -0,0 +1,9 @@ +% etude_fonction_particulier : petit script pour étudier la fonction +% f(x)=2/sqrt(3)*cos(x)+3-2sqrt(3) + +g=inline('2/sqrt(3)*cos(x)+3-2*sqrt(3)'); +gp=inline('-2/sqrt(3)*sin(x)'); +x=0:0.01:pi/3; +plot(x,g(x)-x,x,gp(x)); +legend('courbe de g(x)-x','courbe de g''(x)'); + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4e/sqrtm_iteration_newton.m b/tel/TPmatlab/equation_nonlineaire/TP4e/sqrtm_iteration_newton.m new file mode 100644 index 0000000..47eb3fe --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4e/sqrtm_iteration_newton.m @@ -0,0 +1,70 @@ +function [er,X]=sqrtm_iteration_newton(A,X0,epsilon) + +% +% sqrtm_iteration_newton : calcul d'une racine carrée de matrice par la méthode de Newton. +% +% +% ********************************************************* +% +% [er,X]=sqrtm_iteration_newton(A,X0,espilon) +% calcule les différents itérés de +% la méthode de Newton pour la recherche d'une racine de matrice. +% A chaque itération, elle affiche et stocke le résultat X_k ainsi que l'erreur +% er(k)=|||X_k^2-A|||. +% +% variables d'entrées : +% * A :matrice dont on veut une racine. +% * X0 : matrice initiale +% * epsilon : paramètre de précision +% variables de sortie : +% * er: suite des erreurs, définies par er(k)=norm(X_k^2-A). +% * X : suite des itérés (matriciels) calculés. +% +% ATTENTION : nombre maximum d'itération défini par nmax +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +% nombre maximum d'itération défini par nmax +nmax=100; + +% Contrôles d'entrée +[n,p]=size(A); +if (n~=p) + error('attention, A n''est pas carrée : arrêt du programme'); +end +if (norm(A*X0-X0*A,inf)~=0) + disp('attention, la matrice initiale et A ne commutent pas !!'); +end + + + +% Corps d'algorithme +er=[]; +x=X0; +X=x; +k=0; +test=(norm(X0^2-A,inf)~=0); +while(test) + y=(x+A/x)/2; + disp(['itération numéro ',int2str(k)]); + y2=y^2; + erloc=norm(y2-A); + er=[er erloc]; + X=[X y]; + test=~((norm(y2-A,inf)==0)|(k>=nmax)|(norm(x-y,inf)<=epsilon)); + disp('itéré : '); + disp(y); + disp(['erreur : ',num2str(erloc)]); + disp('appuyez sur une touche pour continuer'); + pause; + k=k+1; + x=y; +end +if (n>=nmax) + disp('attention, n>nmax'); +end diff --git a/tel/TPmatlab/equation_nonlineaire/TP4f/iteration_ordre3_gen.m b/tel/TPmatlab/equation_nonlineaire/TP4f/iteration_ordre3_gen.m new file mode 100644 index 0000000..87464fe --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4f/iteration_ordre3_gen.m @@ -0,0 +1,93 @@ +function [n,X]=iteration_ordre3_gen(x0,nmax,epsilon,fcn) + + +% iteration_ordre3_gen : calcul d'une racine la méthode d'ordre 3 générale +% +% ********************************************************* +% +% [n,X]=iteration_ordre3_gen(x0,nmax,epsilon,fcn) renvoie +% l'approximation de la racine d'une fonction par la méthode d'ordre trois +% générale (les dérivées f' et f'' sont calculées symboliquement) +% +% variables d'entrées : +% * nmax : nombre maximal d'itérations +% * x0 : le premier terme +% * epsilon : précision souhaitée +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% +% variables de sortie : +% * X est une vecteur qui contient les valeurs x_0, x_1, ...,x_n +% * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|=nmax +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% Contrôles d'entrée +qs=exist(fcn); +if (qs~=2) & (qs~=5) & (qs~=1) + error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)'); +end + +% corps d'algorithme +syms vx +g=feval(fcn,vx); +fcnpsymb=diff(g,vx); +fcnppsymb=diff(fcnpsymb,vx); +n=0; +test=(feval(fcn,x0)~=0); +x=x0; +X=x; +fx=feval(fcn,x); +fxp=subs(fcnpsymb,vx,x); +if ~(isnumeric(fxp)); + xpp=eval(fxp); +end +fxpp=subs(fcnppsymb,vx,x); +if ~(isnumeric(fxpp)); + xpp=eval(fxpp); +end +while(test) + A=fxpp/2; + B=fxp; + C=fx; + if (C==0) + d=0; + else + r=sqrt(B^2-4*A*C); + Bp=B+r; + Bm=B-r; + Bpn=abs(Bp); + Bmn=abs(Bm); + if (Bpn=nmax)|(abs(x-xnp1)<=epsilon)); + x=xnp1; + fx=feval(fcn,x); + fxp=subs(fcnpsymb,vx,x); + if ~(isnumeric(fxp)); + xpp=eval(fxp); + end + fxpp=subs(fcnppsymb,vx,x); + if ~(isnumeric(fxpp)); + xpp=eval(fxpp); + end + n=n+1; + X=[X x]; +end +if (n>=nmax) + disp('attention, n>nmax'); +end diff --git a/tel/TPmatlab/equation_nonlineaire/TP4f/iteration_ordre3_gen_symb.m b/tel/TPmatlab/equation_nonlineaire/TP4f/iteration_ordre3_gen_symb.m new file mode 100644 index 0000000..c86ce6d --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4f/iteration_ordre3_gen_symb.m @@ -0,0 +1,102 @@ +function [n,X]=iteration_ordre3_gen_symb(x0,nmax,epsilon,fcn) + + +% iteration_ordre3_gen_symb : calcul d'une racine la méthode d'ordre 3 générale (symbolyque) +% +% ********************************************************* +% +% [n,X]=iteration_ordre3_gen_symb(x0,nmax,epsilon,fcn) renvoie +% l'approximation de la racine d'une fonction par la méthode d'ordre trois +% générale (les dérivées f' et f'' sont calculées symboliquement) +% calcul en symbolique (x0 et fcn peuvent l'être) +% +% variables d'entrées : +% * nmax : nombre maximal d'itérations +% * x0 : le premier terme (éventuellement symbolique) +% * epsilon : précision souhaitée +% * fcn est une chaîne de caractère représentant la fonction (éventuellement symbolique) +% (de type inline, builtin ou par fichier M-file); +% +% variables de sortie : +% * X est une vecteur qui contient les valeurs x_0, x_1, ...,x_n +% * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|=nmax +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% Contrôles d'entrée +qs=exist(fcn); +if (qs~=2) & (qs~=5) & (qs~=1) + error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)'); +end + +% corps d'algorithme +syms vx +g=feval(fcn,vx); +fcnpsymb=diff(g,vx); +fcnppsymb=diff(fcnpsymb,vx); +n=0; +test=(feval(fcn,x0)~=0); +x=x0; +X=x; +fx=feval(fcn,x); +fxp=subs(fcnpsymb,vx,x); +if ~(isnumeric(fxp)); + xpp=eval(fxp); +end +fxpp=subs(fcnppsymb,vx,x); +if ~(isnumeric(fxpp)); + xpp=eval(fxpp); +end +while(test) + A=fxpp/2; + B=fxp; + C=fx; + if (C==0) + d=0; + else + r=sqrt(B^2-4*A*C); + Bp=B+r; + Bm=B-r; + Bpn=abs(Bp); + Bmn=abs(Bm); + if (isnumeric(Bpn)) + if (Bpn=nmax)|(abs(double(x-xnp1))<=epsilon)); + x=xnp1; + fx=feval(fcn,x); + fxp=subs(fcnpsymb,vx,x); + if ~(isnumeric(fxp)); + xpp=eval(fxp); + end + fxpp=subs(fcnppsymb,vx,x); + if ~(isnumeric(fxpp)); + xpp=eval(fxpp); + end + n=n+1; + X=[X x]; +end +if (n>=nmax) + disp('attention, n>nmax'); +end diff --git a/tel/TPmatlab/equation_nonlineaire/TP4g/aitken.m b/tel/TPmatlab/equation_nonlineaire/TP4g/aitken.m new file mode 100644 index 0000000..24de69e --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4g/aitken.m @@ -0,0 +1,124 @@ +function [ait,deb_acc]=aitken(x0,exp_g,eps_seuil,n_max,nul_seuil) +% produit les accélérés d'Aitken d'une suite du point fixe et l'indice de début d'accélération. +% La convergence de la suite initiale est supposée connue, +% ainsi que celle des accélérés. +% +% variables d'entrée +% x0 est la valeur initiale, +% exp_g le nom de la fonction d'itération passé sous la forme 'g(x)' +% eps_seuil désigne comme annoncé dans l'ouvrage le seuil de variation +% relative des itérés qui déclenche le processus d'accélération. +% n_max est la taille maximale du vecteur attendu en sortie; +% on se limite à des valeurs d'itérés distinctes au sens de nul_seuil. +% nul_seuil est un réel positif en-dessous duquel un réel est considéré +% nul; +% plus précisément pour un vecteur d'itérés on ne prend plus en compte +% les valeurs distantes de la valeur finale de moins de nul_seuil. +% +% variables de sortie +% ait désigne le vecteur constitué des premiers itérés ordinaires de x0 +% par la fonction g, +% suivis dès que possible au sens des paramètres passés, +% par les accélérés d'Aitken, +% avec arret dès que les suivants sont identiques au dernier terme +% au sens de nul_seuil choisi. +% deb_acc est l'indice au sens matlab du premier accéléré produit. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions appelées * +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% spécifications +format long e; + + +% tests d'entrée +if nargin~=5 + error('Nombre des arguments incorrect'); +end +if n_max<2 + error('taille n_max insuffisante pour accélérer'); +end +% Il convient sans doute de tester d'autres paramètres, leur type etc... +% selon la dangerosité des utilisateurs! + + +% Algorithme proprement dit + +% Génération des itérés simples par la fonction g; +% on peut remplacer par une écriture matricielle ce bloc. Il sera moins +% lisible. +x=x0;val(1)=x; +for i=1:n_max-1 + x=eval(exp_g); + val(i+1)=x; +end + + +% Mise à l'abri de misères! +% On évite des valeurs identiques à partir d'un certain rang. +ind=find(abs(diff(val))>nul_seuil);val=val(ind); +n_max=size(val,2); + +% recherche des ratios de variation relative des éléments de val; +% voir cours. + +% Ici q(n) représentera ce qui dans le cours est q d'indice n +% contrairement à l'ordinaire matlab. +q=diff(val(1:n_max-1))./diff(val(2:n_max)); + +% variation des qn +rel=diff(q)./q(2:n_max-2); +ind=find(abs(rel)=nul_seuil); +ait=ait(ind); + + +% affichage des résultats à la demande de l'utilisateur +% Ce bloc peut etre remplacé par le passage d'un champ complémentaire +% dans la fonction aitken, par exemple dem_affichage, +% qui serait un booléen égal à 1 +% ssi l'utlisateur veut afficher les résultats. +% Ainsi l'utilisateur par une copie de la fenetre de commande +% après exécution récupère et imprime l'appel de fonction +% avec les paramètres passés, suivi des résultats obtenus. + + +ch=input('Voulez-vous afficher les résultats? (Si oui, taper 1)','s'); +if str2num(ch)==1 + % à modifier selon l'impression souhaitée + nb_lig=3; + + s=[]; + for i=1:size(ait,2) + s=[s sprintf('ait%2.0f = %12.10f ',i,ait(i))]; + % pourrait etre remplacé par un calibrage en fonction + % de nul_seuil du type cal=-log10(2*nul_seuil) + % à passer après transformation chaine à la place de 12.10f. + if rem(i,nb_lig)==0|(i==size(ait,2)) + disp(s);sprintf('\n'); + s=[]; + end + end +end + +% fin de fonction aitken diff --git a/tel/TPmatlab/equation_nonlineaire/TP4g/entrelacs_conv_anim.m b/tel/TPmatlab/equation_nonlineaire/TP4g/entrelacs_conv_anim.m new file mode 100644 index 0000000..2e5c7ec --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4g/entrelacs_conv_anim.m @@ -0,0 +1,156 @@ +function entrelacs_conv_anim(f1,f2) +% Propose une visualisation animée de la convergence d'une suite. +% +% variables d'entrée +% f1 désigne un vecteur de n1 valeurs d'une suite convergente +% dont la dernière est considérée comme la limite. +% L'écart à la limite sera représenté par des rayons polaires +% répartis sur un cercle proportionnellement aux gains de variation relative, +% rapportés à la somme de ceux-ci, tous considérés en valeur absolue. +% Le rayon est l'écart de chaque valeur à la "limite". +% Idem en f2. + + +% Tests d'entrée +if nargin~=2 + error('Nombre des arguments incorrect'); +end + +% ménage! +clf; + +% Calibrage des données: ces constantes modifiables par l'utilisateur +% pourraient aussi etre passées comme champs complémentaires. +% nul_seuil: un nombre inférieur en valeur absolue à nul_seuil est considéré nul. +% rap_seuil: pourcentage de la fenetre en dessous duquel on dilate l'image. +% u: unité de pause en secondes; règle la vitesse de production des images. + + +nul_seuil=0.5E-11; +rap_seuil=0.1; +u=1; + +% détermination des rayons vecteurs, angles +% et initialisations diverses relatives à f1 +r1=abs(f1-f1(size(f1,2)));ind1=find(r1>nul_seuil); +r1=r1(ind1); +nb_pas1=size(r1,2); +att1=zeros(1,nb_pas1);alph1=zeros(1,nb_pas1); +var_rel1=abs(diff(r1))./abs(r1(2:nb_pas1)); +att1(2:nb_pas1)=1./var_rel1*u; +alph1(2:nb_pas1)=2*pi/sum(var_rel1)*cumsum(var_rel1); + +% détermination des rayons vecteurs, angles +% et initialisations diverses relatives à f2 +r2=abs(f2-f2(size(f2,2)));ind2=find(r2>nul_seuil); +r2=r2(ind2); +nb_pas2=size(r2,2); +att2=zeros(1,nb_pas2);alph2=zeros(1,nb_pas2); +var_rel2=abs(diff(r2))./abs(r2(2:nb_pas2)); +att2(2:nb_pas2)=1./var_rel2*u; +alph2(2:nb_pas2)=2*pi/sum(var_rel2)*cumsum(var_rel2); + +% création de la table globale +% qui sera ordonnée temporellement ultérieurement. +date1=cumsum(att1);num1=ones(1,nb_pas1); +X1=r1.*cos(alph1);Y1=r1.*sin(alph1); +date2=cumsum(att2);num2=3*ones(1,nb_pas2); +X2=r2.*cos(alph2);Y2=r2.*sin(alph2); + +date=[date1 date2];num=[num1 num2]; +X=[X1 X2];Y=[Y1 Y2];r=[r1 r2]; +[tdate,ind]=sort(date); +tnum=num(ind); tr=r(ind); +tX=X(ind);tY=Y(ind);tatt=[0 diff(tdate)]; + + +% création des valeurs qui s'enroulent autour de +% la limite avec ajustement d'échelle automatique +% initialisations +coef=1.1*min(r1(1),r2(1)); +lim1=[-coef coef -coef coef]; +c=rap_seuil*coef; +x=[-c,-c,c,c,-c];y=[-c,c,c,-c,-c]; +tot=nb_pas1+nb_pas2;compt1=0;compt2=0; +encor1=1;encor2=1; + +n=1; + + +% préparation des fenetres +figure; +subplot(2,2,1);title('Conv de suite 1');hold on; +subplot(2,2,2);title('Var rel cum pour suite 1');hold on; +subplot(2,2,3);title('Conv de suite 2');hold on; +subplot(2,2,4);title('Var rel cum pour suite 2');hold on; + +max1=max(nb_pas1,nb_pas2); +max2=max(sum(var_rel1),sum(var_rel2)); +lim2=1.1*[2 max1 0 max2]; + +subplot(2,2,2); +axis(lim2);plot((2:nb_pas1),cumsum(var_rel1),'k.-'); +hold on; +subplot(2,2,4); +axis(lim2);plot((2:nb_pas2),cumsum(var_rel2),'k.-'); +hold on; + +while n<=tot + if (1/coef*tr(n)>rap_seuil) + % suivi des points traites + ch=num2str(1+(tnum(n)-1)/2); + ch_x=['compt' ch];ch_y=['var_rel' ch]; + eval([ch_x '=' ch_x '+1;']);tx=eval(ch_x); + ty=eval(['sum(' ch_y '(1:tx-1))']); + + % tracés dans les fenetres + subplot(2,2,tnum(n));axis(lim1); + line([0,tX(n)],[0,tY(n)]);hold on; + subplot(2,2,tnum(n)+1);axis(lim2); + plot([tx],[ty],'cs'); + hold on; + + % prise d'une pause (bien méritée!) + pause(tatt(n)); + + % mise à jour de n + n=n+1; + j=find(tnum(n:tot)==1); + if (j==[])&(encor1) + encor1=0; + subplot(2,2,1);axis(lim1); + plot([0],[0],'rp-'); + legend('Limite ''atteinte''',0); + hold on; + end + j=find(tnum(n:tot)==3); + if (j==[])&(encor2); + encor2=0; + subplot(2,2,3);axis(lim1); + plot([0],[0],'rp-'); + legend('Limite ''atteinte''',0); + hold on; + end + + else + % annonce d'agrandissement + c=rap_seuil*coef; + x=[-c,-c,c,c,-c];y=[-c,c,c,-c,-c]; + if encor1 + subplot(2,2,1);axis(lim1); + plot(x,y,'r.--'); + hold on;pause(u/4); + end + if encor2 + subplot(2,2,3);axis(lim1); + plot(x,y,'r.--'); + hold on;pause(u/4); + end + + % changement d'échelle + k=1.1/coef*tr(n); + lim1=k*lim1;coef=k*coef; + end +end + +% fin de fonction \ No newline at end of file diff --git a/tel/TPmatlab/equation_nonlineaire/TP4g/point_fixe.m b/tel/TPmatlab/equation_nonlineaire/TP4g/point_fixe.m new file mode 100644 index 0000000..71a0927 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4g/point_fixe.m @@ -0,0 +1,51 @@ +function val=point_fixe(x0,nom_g,n_max) +% produit un vecteur d'itérés d'une valeur initiale par une fonction g +% +% variable d'entrée +% x0 est la valeur initiale choisie; +% nom_g l'expression de la fonction donnée +% sous la forme d'une chaine 'g(x)'; +% n_max est un entier naturel donnant la longueur du vecteur des itérés. +% +% variables de sortie +% val est le vecteur des itérés de x0 par g de longueur n_max. + +% On pourra tester cette fonction en tapant en fenetre de commande +% point_fixe(2,'(60*x-36-6*x^2)^(1/4)',20). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions appelées * +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% spécifications +format long e; + + +x=x0;val(1)=x; +for i=1:n_max-1 + x=eval(nom_g); + val(i+1)=x;val; +end + +ch=input('Voulez-vous afficher les résultats? (Si oui, taper 1)','s'); +if str2num(ch)==1 + % à modifier selon l'impression souhaitée + nb_lig=3; + + s=[]; + for i=1:size(val,2) + s=[s sprintf('ptf%2.0f = %12.10f ',i,val(i))]; + % pourrait etre remplacé par un calibrage en fonction + % de nul_seuil du type cal=-log10(2*nul_seuil) + % à passer après transformation chaine à la place de 12.10f. + if (rem(i,nb_lig)==0)|(i==size(val,2)) + disp(s);sprintf('\n'); + s=[]; + end + end +end + +% fin de point_fixe + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4g/steffensen.m b/tel/TPmatlab/equation_nonlineaire/TP4g/steffensen.m new file mode 100644 index 0000000..268d3cc --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4g/steffensen.m @@ -0,0 +1,178 @@ +function [stef,deb_acc]=steffensen(x0,exp_g,eps_seuil,n_max,nul_seuil) +% Produit accélérés de Steffensen d'une suite du point fixe et indice de +% début d'accélération. +% La convergence de la suite initiale est supposée connue, +% ainsi que celle des accélérés. +% Ce fichier bénéficie de la connaissance préalable de la fonction ait. +% +% variables d'entrée +% x0 est la valeur initiale, +% exp_g le nom de la fonction d'itération passé sous la forme 'g(x)' +% eps_seuil désigne comme annoncé dans l'ouvrage le seuil de variation +% relative des itérés, +% qui déclenche le processus d'accélération. +% n_max est la taille maximale du vecteur attendu en sortie; +% on se limite à des valeurs d'itérés distinctes au sens de nul_seuil. +% nul_seuil est un réel positif en-dessous duquel +% un réel est considéré nul. +% +% variables de sortie +% stef désigne le vecteur constitué des premiers itérés ordinaires de x0 +% par la fonction g, +% suivis dès que possible au sens des paramètres passés, +% par les accélérés de Steffensen, +% avec arret dès que les suivants sont identiques au sens de nul_seuil +% choisi. +% deb_acc est l'indice au sens matlab du premier accéléré produit. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions appelées * +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% spécifications +format long e; + + +% Tests d'entrée +if nargin~=5 + error('Nombre des arguments incorrect'); +end +if n_max<2 + error('taille n_max insuffisante pour accélérer'); +end +% Il convient sans doute de tester d'autres paramètres, leur type etc... +% selon la dangerosité des utilisateurs! +% Test complémentaire sur donnée +if n_max<2 + error('taille n_max insuffisante pour accélérer'); +end + +% Algorithme proprement dit + +% Initialisation des vecteurs stef et util + +stef=zeros(1,n_max);stef(1)=x0; +% util vecteur de longueur 4, contient initialement les itérés lents +% ordinaires. +x=x0;util(1)=x; +for i=1:3 + x=eval(exp_g); + util(i+1)=x; +end + +% Initialisations +% numéro du premier candidat à l'accélération +% au sens des indices matlab +compt=2;deb_acc=2; + +% Initialisations avant l'entrée en boucle +% en pratique toujours parcourue une fois... + +q_util=diff(util(1:3))./diff(util(2:4)); +rel=diff(q_util)/q_util(2); +% définition du booléen d'autorisation d'accélération: voir cours. +b_acc=(abs(rel)nul_seuil) + + % calcul de l'accéléré suivant + % ici par formule du delta2 d'Aitken! + val=acc(3)-(diff(acc(2:3))^2)/diff(diff(acc)); + stef(compt)=val;acc(3)=val; + + % Il convient de noter que les accélérés obtenus pourraient etre + % améliorés en considérant les deux accélérés obtenus à partir de + % trois itérés successifs; ceci éviterait le phénomène qui apparait + % sous la fonction entrelacs... appliquée aux steffensen. + + % préparation du tour suivant + acc(1:2)=acc(2:3); + x=acc(2);acc(3)=eval(exp_g); + compt=compt+1; +end + +% suppression des valeurs rendues inutiles après l'arret du calcul +stef;ind=find(stef~=0);stef=stef(ind); + + +% affichage des résultats à la demande de l'utilisateur +% Ce bloc peut etre remplacé par le passage d'un champ complémentaire +%dans la fonction steffensen par exemple dem_affichage, +% qui serait un booléen égal à 1 +% ssi l'utilisateur veut afficher les résultats. +% Ainsi l'utilisateur par une copie de la fenetre de commande +% après exécution récupère +% et imprime l'appel de fonction avec les paramètres passés +% suivi des résultats obtenus. + +ch=input('Voulez-vous afficher les résultats? (Si oui, taper 1)','s'); +if str2num(ch)==1 + % à modifier selon l'impression souhaitée + nb_lig=3; + + s=[]; + for i=1:size(stef,2) + s=[s sprintf('stef%2.0f = %12.10f ',i,stef(i))]; + % pourrait etre remplacé par un calibrage en fonction + % de nul_seuil du type cal=-log10(2*nul_seuil) + % à passer après transformation chaine à la place de 12.10f. + if rem(i,nb_lig)==0|(i==size(stef,2)) + disp(s);sprintf('\n'); + s=[]; + end + end +end + +% fin de fonction steffensen + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4g/visu1_conv_anim.m b/tel/TPmatlab/equation_nonlineaire/TP4g/visu1_conv_anim.m new file mode 100644 index 0000000..fbca55e --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4g/visu1_conv_anim.m @@ -0,0 +1,115 @@ +function visu1_conv_anim(f1) +% Propose une visualisation animée de la convergence d'une suite. +% +% variables d'entrée +% f1 désigne un vecteur de n1 valeurs d'une suite convergente +% dont la dernière est considérée comme la limite. +% L'écart à la limite sera représenté par des rayons polaires +% répartis sur un cercle proportionnellement aux gains de variation relative, +% rapportés à la somme de ceux-ci, tous considérés en valeur absolue. +% Le rayon est l'écart de chaque valeur à la "limite". + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions appelées * +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Tests d'entrée +if nargin~=1 + error('Nombre des arguments incorrect'); +end + +% ménage! +clf; + +% Calibrage des données: ces constantes modifiables par l'utilisateur +% pourraient aussi etre passées comme champs complémentaires. + +% nul_seuil: un nombre inférieur en valeur absolue à nul_seuil est considéré nul. +% rap_seuil: pourcentage de la fenetre en dessous duquel on dilate l'image. +% u unité de temps en seconde, pour l'expression d'une pause +% déterminant la vitesse de défilement des images. +nul_seuil=0.5E-11; +rap_seuil=0.1; +u=2; + + +% détermination des rayons vecteurs, angles +% et initialisations diverses relatives à f1 + +r1=abs(f1-f1(size(f1,2)));ind1=find(r1>nul_seuil); +r1=r1(ind1); +nb_pas1=size(r1,2); +att1=zeros(1,nb_pas1);alph1=zeros(1,nb_pas1); +var_rel=abs(diff(r1))./abs(r1(2:nb_pas1)); +att1(2:nb_pas1)=1./var_rel*u; +alph1(2:nb_pas1)=2*pi/sum(var_rel)*cumsum(var_rel); +coef=1.1*r1(1); +lim=[-coef coef -coef coef]; + + +% création des valeurs qui s'enroulent autour de +% la limite avec ajustement d'échelle automatique + +n=1; + +% préparation de la fenetre +title('Visualisation de convergence'); +c=rap_seuil*coef; +X=[-c,-c,c,c,-c];Y=[-c,c,c,-c,-c]; + +while n<=nb_pas1 + if (1/coef*r1(n)>rap_seuil) + + x1=r1(n)*cos(alph1(n)); + y1=r1(n)*sin(alph1(n)); + % tracé + axis(lim); + line([0,x1],[0,y1]);hold on; + % prise d'une pause (bien méritée!) + pause(att1(n)); + % mise à jour de n + n=n+1; + else + % annonce d'agrandissement + c=rap_seuil*coef; + X=[-c,-c,c,c,-c];Y=[-c,c,c,-c,-c]; + axis(lim); + plot(X,Y,'r.--'); + hold on;pause(u/4) + % changement d'échelle + k=1.1/coef*r1(n); + lim=k*lim;coef=k*coef; + end +end + +% indication de fin de convergence +axis(lim); +plot([0],[0],'rp-');legend('Limite ''atteinte''',0); +hold off; + +% fin de fonction + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4g/visu2_conv_anim.m b/tel/TPmatlab/equation_nonlineaire/TP4g/visu2_conv_anim.m new file mode 100644 index 0000000..941c575 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4g/visu2_conv_anim.m @@ -0,0 +1,154 @@ +function visu2_conv_anim(f1,f2) +% Propose une visualisation animée de la convergence comparée de deux suites. +% +% variables d'entrée +% f1 désigne un vecteur de n1 valeurs qui seront représentées en polaire +% par pas angulaire de 2*pi/n1; +% le rayon est l'écart de chaque valeur à la dernière, considérée comme limite. +% Idem en f2. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions appelées * +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% Tests d'entrée +if nargin~=2 + error('Nombre des arguments incorrect'); +end + +% ménage! +clf; + +% Calibrage des données: ces constantes modifiables par l'utilisateur +% pourraient aussi etre passées comme champs complémentaires. +% nul_seuil: un nombre inférieur en valeur absolue à nul_seuil est considéré nul. +% rap_seuil: pourcentage de la fenetre en dessous duquel on dilate l'image. +% u: unité de pause en secondes; règle la vitesse de production des images. + + +nul_seuil=0.5E-11; +rap_seuil=0.1; +u=1; + +% détermination des rayons vecteurs +r1=abs(f1-f1(size(f1,2)));ind1=find(r1>nul_seuil); +r1=r1(ind1); +r2=abs(f2-f2(size(f2,2)));ind2=find(r2>nul_seuil); +r2=r2(ind2); + +% calibrage des axes et choix du pas angulaire +alph1=2*pi/size(r1,2); +alph2=2*pi/size(r2,2); +coef=1.1*min(r1(1),r2(1)); +lim=[-coef coef -coef coef]; + +% création des valeurs qui s'enroulent autour de +% la limite avec ajustement d'échelle automatique + +tailles=[size(r1,2),size(r2,2)]; +nb_pas1=min(tailles);nb_pas2=max(tailles); +n=1; + +% préparation des fenetres +subplot(2,1,1);title('Convergence de suite 1'); +subplot(2,1,2);title('Convergence de suite 2'); +% on peut faire plus rusé pour le titre! +c=rap_seuil*coef; +X=[-c,-c,c,c,-c];Y=[-c,c,c,-c,-c]; + +while n<=nb_pas1 + if (1/coef*min(r1(n),r2(n))>rap_seuil) + + x1=r1(n)*cos((n-1)*alph1); + y1=r1(n)*sin((n-1)*alph1); + x2=r2(n)*cos((n-1)*alph2); + y2=r2(n)*sin((n-1)*alph2); + + % tracé dans les deux fenetres + subplot(2,1,1);axis(lim); + line([0,x1],[0,y1]);hold on; + subplot(2,1,2);axis(lim); + line([0,x2],[0,y2]);hold on; + + % prise d'une pause (bien méritée!) + pause(2*u); + + % mise à jour de n + n=n+1; + else + % annonce d'agrandissement + c=rap_seuil*coef; + X=[-c,-c,c,c,-c];Y=[-c,c,c,-c,-c]; + subplot(2,1,1);axis(lim); + plot(X,Y,'r.--'); + hold on;pause(2) + subplot(2,1,2);axis(lim); + plot(X,Y,'r.--'); + hold on;pause(2); + % changement d'échelle + k=1.1/coef*min(r1(n),r2(n)); + lim=k*lim;coef=k*coef; + end +end +% Variante : il est possible d'utiliser +% la fonction movie et de créer une animation + +% indication de fin de convergence + +if nb_pas1~=nb_pas2 + ind=find(tailles==nb_pas2); + r3=eval(['r' num2str(ind)]); + alph3=eval(['alph' num2str(ind)]); + num=3-ind; + subplot(2,1,num);axis(lim); + plot([0],[0],'rp-');legend('Limite ''atteinte''',0); + hold on; + % préparation de fin de tracé + num=3-num; +else + subplot(2,1,1);axis(lim); + plot([0],[0],'rp-');legend('Limite ''atteinte''',0); + hold on;subplot(2,1,2);axis(lim); + plot([0],[0],'bp-');legend('Limite ''atteinte''',0); + hold on; +end + +% tracé de la courbe restante éventuelle + +while n<=nb_pas2 + if (1/coef*r3(n)>rap_seuil) + + x3=r3(n)*cos((n-1)*alph3); + y3=r3(n)*sin((n-1)*alph3); + + % tracé dans la fenetre convenable + subplot(2,1,num);axis(lim); + line([0,x3],[0,y3]);hold on; + + % prise d'une pause (bien méritée!) + pause(2); + + % mise à jour de n + n=n+1; + else + % annonce d'agrandissement + c=rap_seuil*coef; + X=[-c,-c,c,c,-c];Y=[-c,c,c,-c,-c]; + subplot(2,1,num);axis(lim); + plot(X,Y,'r.--'); + hold on;pause(2); + % changement d'échelle + k=1.1/coef*r3(n); + lim=k*lim;coef=k*coef; + end +end + +% limite atteinte +subplot(2,1,num);axis(lim); +plot([0],[0],'bp-');legend('Limite ''atteinte''',0); +hold off; + +% fin de fonction \ No newline at end of file diff --git a/tel/TPmatlab/equation_nonlineaire/TP4h/affiche_racines.m b/tel/TPmatlab/equation_nonlineaire/TP4h/affiche_racines.m new file mode 100644 index 0000000..1d96b8d --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4h/affiche_racines.m @@ -0,0 +1,33 @@ +function affiche_racines(ens_rac,exp_pol3,f_seuil) +% affiche les racines sous différentes formes possibles +% +% variables d'entrée +% ens_rac est le vecteur des racines de l'équation traitée; +% exp_pol3 est la chaine associée à pol3; +% f_seuil désigne le module de pol3(x) pour lequel il est assimilé à zéro. +% +% variables de sortie +% on peut extraire les objets fabriqués. + +% tests à écrire; pas de fonctions sous-jacente... + +% au fait +% sortie des résultats sous forme symbolique +disp('Sortie des résultats sous forme symbolique proposée par matlab'); +symb=sym(ens_rac,'r');disp('ens_rac = ');disp (symb); + + +% tentative d'arrondi des racines au sens de f_seuil +% pourrait s'écrire matriciellement ou etre une routine annexe déclenchée +% par le passage d'un champ complémentaire. + +for k=1:size(ens_rac,2) + x=round(ens_rac(k));val=eval(exp_pol3); + if abs(val)0)) + error('Champ init mal écrit. Voir controle'); +end +if (size(arret,1)~=1)|(size(arret,2)~=3) + error('Champ arret mal écrit. Voir controle'); +end +if (size(pol3,1)~=1)|(isstr(pol3)) + error('Champ pol3 mal écrit. pol3 est un vecteur de réels ici'); +end +% compléments éventuels + + +% Algorithme proprement dit + + +% initialisations et préparation des données +ens_rac=[];nb_iter=[]; + +while pol3(1)==0 + pol3(1)=[]; +end +% si on a passé sciemment le vecteur nul, +% on peut s'interroger sur la nullité de l'utilisateur! +deg=size(pol3,2)-1; + +% création de la chaine d'appel de muller_elem associée à pol3 ou autre! +vect=pol3; +exp_f=vect2str(vect); + + +while deg>1 + + % adjonction des nouvelles racines issues de Muller élémentaire + [rac,iter]=muller_elem(init,exp_f,arret); + ens_rac=[ens_rac rac];nb_iter=[nb_iter iter]; + + % deflation + taille=size(rac,2); + switch taille + case 1 + div=[1 -rac];deg=deg-1; + case 2 + div=[1 -sum(rac) prod(rac)];deg=deg-2; + end + + % préparation de l'itération suivante + vect=deconv(vect,div); + exp_f=vect2str(vect); +end + +% fin de déflation éventuelle +if deg==1 + rac=-vect(2)/vect(1); + iter=0; + ens_rac=[ens_rac rac];nb_iter=[nb_iter iter]; +end + +% affichage éventuel des résultats +% à libérer du commentaire éventuellement +% affiche_racines(ens_rac,vect2str(pol3),arret(1)); +% fin de fonction \ No newline at end of file diff --git a/tel/TPmatlab/equation_nonlineaire/TP4h/muller_elem.m b/tel/TPmatlab/equation_nonlineaire/TP4h/muller_elem.m new file mode 100644 index 0000000..5d30ada --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4h/muller_elem.m @@ -0,0 +1,114 @@ +function [rac,nb_iter]=muller_elem(init,exp_f,arret) +% Produit une racine d'équation (E) f(x)=0 par méthode de Muller. +% Le calcul se termine par un forçage à une écriture symbolique de rac, +% qui est affichée en plus de la numérique! +% Ceci peut etre aisément modifié! +% +% variables d'entrée +% init est un vecteur de trois valeurs initiales distinctes +% pour la méthode de Muller classique. On peut généraliser! +% exp_f est l'expression de la fonction f sous la forme 'f(x)'; +% arret est un vecteur de trois réels [f_seuil,nul_seuil,n_max] +% dont la signification est la suivante: +% si |f(x)|0)) + error('Champ init mal écrit. Voir controle'); +end +if (size(arret,1)~=1)|(size(arret,2)~=3) + error('Champ arret mal écrit. Voir controle'); +end +% compléments éventuels + + +% Algorithme proprement dit + +% Etape 1:initialisations +X=init;Y=zeros(1,3); + +% La suite serait avantageusement remplacée par une écriture vectorielle +% mais l'utilisateur calculera faux sans le voir..; Il faudrait écrire +% 'x.^2+1' au lieu de 'x^2+1' +for k=1:3 + x=X(k);Y(k)=eval(exp_f); +end + +d=diff(X);dif_div=diff(Y)./d; +nb_iter=0; +bool_arret=0; + +% Etape 2 +while bool_arret==0 + nb_iter=nb_iter+1; + % détermination du polynome d'interpolation de f sur le support x + % et résolution de l'équation approchée + An=diff(dif_div)/sum(d); + Bn=dif_div(2)+An*d(2); + x=X(3);Cn=eval(exp_f); + if abs(Cn)arret(3))|(abs(ajout)=0 + s='+'; +end +exp_f=[s num2str(vect(deg+1))]; +for k=1:deg + s=''; + if vect(deg+1-k)>=0 + s='+'; + end + exp_f=[s num2str(vect(deg+1-k)) '*x^' num2str(k) exp_f]; +end + +% fin de fonction diff --git a/tel/TPmatlab/equation_nonlineaire/TP4i/partie_1/fonct1.m b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_1/fonct1.m new file mode 100644 index 0000000..9018679 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_1/fonct1.m @@ -0,0 +1,9 @@ +function z=fonct1(omega,t) +% calcule pour des vecteurs convenables z=F(omega,t) dans le tp4I. +% variables d'entrée +% omega et t sont définis dans l'énoncé du tp et dans trace_surface +% +% variables de sortie +% z est la valeur définie dans l'énoncé du tp4I. + +z=2*(t.*omega-pi).*sin(omega.*t)+cos(omega.*t); diff --git a/tel/TPmatlab/equation_nonlineaire/TP4i/partie_1/test_oscillation.m b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_1/test_oscillation.m new file mode 100644 index 0000000..27e349b --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_1/test_oscillation.m @@ -0,0 +1,78 @@ +function [chaine]=test_oscillation(x0) +% calcule les itérés de Newton de la valeur t issue de l'étude de la +% question 1 du tp4I. +% +% variable d'entrée +% x0 représente la solution de l'équation F(1,t)=0 qui devrait fournir +% une suite d'itérés de Newton oscillante. +% +% variable de sortie +% annonce ce qui s'est passé effectivement lors du calcul +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% tests à compléter éventuellement +if nargin~=1 + error('Nombre des arguments incorrect'); +end +clf;format long e; + +% définitions de constantes modifiables par l'utilisateur + +nul_seuil=1E-10; +iter_max=30; + + +% test complémentaire. On peut passer val_omega comme champ pour détourner +% cette fonction de sa seule application au cas omega=1. +val_omega=1; % valeur non nulle, attention! +if (x0<=0)|(x0>=2*pi/val_omega) + error('Le champ x0 est incorrect. Voir test complémentaire'); +end + + +% initialisations +v=zeros(1,iter_max+1); +nb_iter=1; +v(1)=x0;v(2)=v(1)+cos(v(1))/sin(v(1));x1=v(2); +normal=1; + +% boucle générale +while nb_iter=19 + + % L'utilisateur n'a rien choisi! + disp('Trop tard. On sort !'); + end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +sortie; + +% fin de demo! \ No newline at end of file diff --git a/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/det_choix.m b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/det_choix.m new file mode 100644 index 0000000..25af552 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/det_choix.m @@ -0,0 +1,22 @@ +function choix=det_choix(num) +% trouve le choix de l'utilisateur pour passer ce champ à traitement_dem +% +% variable d'entrée +% num est un entier qui désigne le numéro du bouton radio +% de creat_fen2_dem ou 0 si on sort du poussoir du meme fichier. +% on ne controle rien! +global h_f2 BR +global choix rep att2 + + +att2=0; +switch num +case 0, + return; +otherwise + if (get(BR(num),'value')==1) + rep=1;choix=num; + end +end + +% fin de fonction diff --git a/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/etude_hypotheses.m b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/etude_hypotheses.m new file mode 100644 index 0000000..8ef0d41 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/etude_hypotheses.m @@ -0,0 +1,50 @@ +function [ens_H,H2,H3,H4,H5]=etude_hypotheses(A,B) +% etudie la validité de l'hypothese num du théorème pour A et B +% +% variables d'entrée +% A et B désignent respectivement les origine et extrémité de l'intervalle +% [a,b] sur lequel on étudie la convergence. +% H2,H3,H4,H5,ens_H sont les matrices booléennes des indices +% de la matrice des couples (a,b) qui vérifient les hypothèses +% (de numéro indiqué pour les premières, +% de conjonction des précédentes pour ens_H) +% du théorème de condition suffisante de convergence +% de la suite de Newton sur [a,b]. + +% tests sommaires +if nargin~=2 + error('Nombre des arguments incorrect'); +end + +% Donnée interne de pol: attention cette donnée agit fortement +% sur l'écriture des contraintes; on évitera de passer pol comme champ. +pol=[1 0 6 -60 36]; + +% étude hypothèse 1 +% inutile: pol est polynome donc C2 sur tout intervalle de R. + +% étude hypothèse 2 +x=A;f_A=eval(vect2str_mat(pol)); +x=B;f_B=eval(vect2str_mat(pol)); +H2=((A0 ici; +% pour l'étude de la racine proche de 2, ce serait pol'(a)<0. +d1_pol=polyder(pol); +x=A;d1_f_A=eval(vect2str_mat(d1_pol)); +x=B;d1_f_B=eval(vect2str_mat(d1_pol)); +H3=((A0)); + +% étude hypothèse 4 +% La dérivée seconde est certainement positive. +% d2_pol=polyder(d1_pol); +H4=(A0) % donc idem en B! + a_opt=min(tempA); + ind_a=(tempA==a_opt); + b_opt=max(tempB(ind_a)); + plot(a_opt,b_opt,'r*'); + h_leg=legend(['Solution optimale: a= ' ... + num2str(a_opt) '; b= ' num2str(b_opt)],0); + % détails + enf=get(h_leg(1),'Children'); + set(enf(1),'Color','r','Marker','*'); + delete(enf(2)); + + end + end + +end + +% complément dans le cas de l'ensemble des hypothèses + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% possibilité de lecture de coordonnées: supprimer le commentaire. +% on clique gauche pour lire; +% on sort en cliquant droit! + +%lect_courb(1); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ouverture de fenetre de choix 3 +creat_fen3_dem; + +% attente jusqu'a 120 unités de temps: unite_att +compt=0;att3=1; + +while (att3==1)&(compt<120) + pause(unite_att); + compt=compt+1; +end + +% fermeture de fenetre 3, non visibilité et nettoyage de fenetre 1 +close(h_f3);close(h_f1); + + +% gestion de réponse ou de son absence +if compt>=120 + rep=0; + disp('Trop tard. Assez attendu.'); + disp('Vous pouvez modifier le temps sous traitement_dem ligne 123'); +end + + +% fin de fonction + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/vect2str_mat.m b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/vect2str_mat.m new file mode 100644 index 0000000..510f919 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4i/partie_2/vect2str_mat.m @@ -0,0 +1,34 @@ +function [exp_f]=vect2str_mat(vect) +% transforme un vecteur d'écriture de polynome en chaine à évaluer +% sous forme matricielle. +% +% variables d'entrée +% vecteur de définition du polynome considéré +% +% variable de sortie +% chaine d'écriture anx^n+...+a0 du polynome considéré +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% controles inexistants ici car opérés en amont + +deg=size(vect,2)-1; + +% création de la chaine d'appel de muller_elem associée à vect +s=''; +if vect(deg+1)>=0 + s='+'; +end +exp_f=[s num2str(vect(deg+1))]; +for k=1:deg + s=''; + if vect(deg+1-k)>=0 + s='+'; + end + exp_f=[s num2str(vect(deg+1-k)) '.*(x.^' num2str(k) ')' exp_f]; +end + +% fin de fonction diff --git a/tel/TPmatlab/equation_nonlineaire/TP4k/calcul_coeff_ordre4.m b/tel/TPmatlab/equation_nonlineaire/TP4k/calcul_coeff_ordre4.m new file mode 100644 index 0000000..0e915f8 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4k/calcul_coeff_ordre4.m @@ -0,0 +1,42 @@ +function [a,b,c,d]=calcul_coeff_ordre4(A) + +% calcul_coeff_ordre4 : mis au point en symbolique d'une méthode d'ordre 4. +% +% [a,b,c,d]=calcul_coeff_ordre4(A) calcule les coefficents a,b,c,d tels que +% la fonction g(x)=a*x+b/x+c/x^3+d/x^5 définisse une suite du point fixe +% avec un ordre de convergence égal à 4 autour de l=sqrt(A); +% +% paramètres d'entrées : A +% paramètres de sorties : a,b,c,d +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% tests préliminaire +if isnumeric(A) + if (A<=0) + error('attention, l''argument doit être strictement positif'); + end +end +% corps d'algorithme +syms x a b c d; +g=a*x+b/x+c/x^3+d/x^5; +gp=diff(g,x); +gp2=diff(g,x,2); +gp3=diff(g,x,3); +B=solve(g-x,gp,gp2,gp3,a,b,c,d); +a=sym(subs(B.a,x,sqrt(sym(A)))); +b=sym(subs(B.b,x,sqrt(sym(A)))); +c=sym(subs(B.c,x,sqrt(sym(A)))); +d=sym(subs(B.d,x,sqrt(sym(A)))); + + + + diff --git a/tel/TPmatlab/equation_nonlineaire/TP4l/calcul_coeff_racine3.m b/tel/TPmatlab/equation_nonlineaire/TP4l/calcul_coeff_racine3.m new file mode 100644 index 0000000..ac5fd69 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4l/calcul_coeff_racine3.m @@ -0,0 +1,36 @@ +function [a,b,c]=calcul_coeff_racine3(A) + +% calcul_coeff_racine3 : mis au point en symbolique d'une méthode d'ordre 3 convergeant vers l^(1/3). +% +% [a,b,c]=calcul_coeff_racine3(A) calcule les coefficents a,b,c tels que +% la fonction g(x)=a*x+b/x^2+c/x^5 définisse une suite du point fixe +% avec un ordre de convergence égal à 3 autour de l=(A)^(1/3); +% +% paramètres d'entrées : A +% paramètres de sorties : a,b,c +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% tests préliminaire +if isnumeric(A) + if (A<=0) + error('attention, l''argument doit être strictement positif'); + end +end +% corps d'algorithme +syms x a b c; +g=a*x+b/x^2+c/x^5; +gp=diff(g,x); +gp2=diff(g,x,2); +B=solve(g-x,gp,gp2,a,b,c); +a=sym(subs(B.a,x,(sym(A))^(1/3))); +b=sym(subs(B.b,x,(sym(A))^(1/3))); +c=sym(subs(B.c,x,(sym(A))^(1/3))); diff --git a/tel/TPmatlab/equation_nonlineaire/TP4m/sqrtm_iteration_ordre3.m b/tel/TPmatlab/equation_nonlineaire/TP4m/sqrtm_iteration_ordre3.m new file mode 100644 index 0000000..ae5c451 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4m/sqrtm_iteration_ordre3.m @@ -0,0 +1,72 @@ +function [er,X]=sqrtm_iteration_ordre3(A,X0,epsilon) + +% +% sqrtm_iteration_ordre3 : calcul d'une racine carrée de matrice par la méthode de Newton. +% +% +% ********************************************************* +% +% [er,X]=sqrtm_iteration_ordre3(A,X0,epsilon) +% calcule les différents itérés de +% la méthode d'ordre 3 pour la recherche d'une racine de matrice. +% A chaque itération, elle affiche et stocke le résultat X_k ainsi que l'erreur +% er(k)=|||X_k^2-A|||. +% +% variables d'entrées : +% * A :matrice dont on veut une racine. +% * X0 : matrice initiale +% * epsilon : paramètre de précision +% variables de sortie : +% * er: suite des erreurs, définies par er(k)=norm(X_k^2-A). +% * X : suite des itérés (matriciels) calculés. +% +% ATTENTION : nombre maximum d'itération défini par nmax +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +% nombre maximum d'itération défini par nmax +nmax=100; + +% Contrôles d'entrée +[n,p]=size(A); +if (n~=p) + error('attention, A n''est pas carrée : arrêt du programme'); +end +sqrtmA=sqrtm(A); +if (norm(sqrtmA*X0-X0*sqrtmA,inf)~=0) + disp('attention, la matrice initiale et la racine de A ne commutent pas !!'); +end + + +% Corps d'algorithme +er=[]; +X=X0; +x=X0; +k=0; +A2=A^2; +test=(norm(x^2-A,inf)~=0); +while(test) + z=inv(x); + y=(3/8)*x+(3/4)*A*z-(1/8)*A2*(z^3); + disp(['itération numéro ',int2str(k)]); + y2=y^2; + erloc=norm(y2-A); + er=[er erloc]; + X=[X y]; + test=~((norm(y2-A,inf)==0)|(k>=nmax)|(norm(x-y,inf)<=epsilon)); + disp('itéré : '); + disp(y); + disp(['erreur : ',num2str(erloc)]); + disp('appuyez sur une touche pour continuer'); + pause; + k=k+1; + x=y; +end +if (n>=nmax) + disp('attention, n>nmax'); +end diff --git a/tel/TPmatlab/equation_nonlineaire/TP4n/pic.m b/tel/TPmatlab/equation_nonlineaire/TP4n/pic.m new file mode 100644 index 0000000..6e069b9 --- /dev/null +++ b/tel/TPmatlab/equation_nonlineaire/TP4n/pic.m @@ -0,0 +1,10 @@ +function res=pic(x) + +% pic(x) est la fonction égale à +% 2*x sur [0,1/2] ; +% 2*(1-x) sur [1/2,1]. + +res=2*x.*(double(x)<=1/2)+2*(1-x).*(double(x)>1/2); +if ~(isnumeric(x)) + res=simplify(res); +end \ No newline at end of file diff --git a/tel/TPmatlab/erreur_algorithmique/TP1d/sommeA.m b/tel/TPmatlab/erreur_algorithmique/TP1d/sommeA.m new file mode 100644 index 0000000..9193a08 --- /dev/null +++ b/tel/TPmatlab/erreur_algorithmique/TP1d/sommeA.m @@ -0,0 +1,40 @@ +function res=sommeA(p,n) + +% sommeA : calcul de somme (version A). +% +% ********************************************************* +% +% res=sommeA(p,n) +% +% calcul de la somme \sum_{i=1}^n 1/(i^p) en utilisant une boucle +% et en sommant dans l'ordre des indices croissants. +% +% variables d'entrées : +% * p,n : deux entiers +% +% variables de sortie +% * res: résultats de la somme +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels + + +% Corps d'algorithme +s=1; +for i=2:n + s=s+1/(i^p); +end +res=s; diff --git a/tel/TPmatlab/erreur_algorithmique/TP1d/sommeB.m b/tel/TPmatlab/erreur_algorithmique/TP1d/sommeB.m new file mode 100644 index 0000000..ebb2c4b --- /dev/null +++ b/tel/TPmatlab/erreur_algorithmique/TP1d/sommeB.m @@ -0,0 +1,36 @@ +function res=sommeB(p,n) + +% sommeB : calcul de somme (version B). +% +% ********************************************************* +% +% res=sommeB(p,n) +% +% calcul de la somme \sum_{i=1}^n 1/(i^p) en utilisant un calcul +% matriciel et en sommant dans l'ordre des indices croissants. +% +% variables d'entrées : +% * p,n : deux entiers +% +% variables de sortie +% * res: résultats de la somme +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels + + +% Corps d'algorithme +res=sum(1./((1:n).^p)); diff --git a/tel/TPmatlab/erreur_algorithmique/TP1d/sommeC.m b/tel/TPmatlab/erreur_algorithmique/TP1d/sommeC.m new file mode 100644 index 0000000..5772c87 --- /dev/null +++ b/tel/TPmatlab/erreur_algorithmique/TP1d/sommeC.m @@ -0,0 +1,36 @@ +function res=sommeC(p,n) + +% sommeC : calcul de somme (version C). +% +% ********************************************************* +% +% res=sommeC(p,n) +% +% calcul de la somme \sum_{i=1}^n 1/(i^p) en utilisant un calcul +% matriciel et en sommant dans l'ordre des indices décroissants. +% +% variables d'entrées : +% * p,n : deux entiers +% +% variables de sortie +% * res: résultats de la somme +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels + + +% Corps d'algorithme +res=sum(1./((n:-1:1).^p)); diff --git a/tel/TPmatlab/erreur_algorithmique/TP1e/affiche_erreur_derivee.m b/tel/TPmatlab/erreur_algorithmique/TP1e/affiche_erreur_derivee.m new file mode 100644 index 0000000..90ab623 --- /dev/null +++ b/tel/TPmatlab/erreur_algorithmique/TP1e/affiche_erreur_derivee.m @@ -0,0 +1,52 @@ +function affiche_erreur_derivee(hm,hM,R) +% +% +% affiche_erreur_derivee : affichage de l'erreur (logarithmique) pour calcul numérique de dérivée. +% +% ********************************************************* +% +% affiche_erreur_derivee(hm,hM,R) : +% affiche le nuage de point (log10(h_i),log10(epsilon(h_i))) +% où pour tout entier i appartenant à (0,...,R), +% h_i=h=hm.^((R-i)/R).*hM.^(i/R); +% epsilon(h_i) est l'erreur commise en remplaçant la dérivée +% de la fonction exponentielle en 1 par (exp(1+h_i)-exp(t-h_i))/(2h_i) +% +% +% variables d'entrées : +% * hm et hM sont deux pas strictement positifs (hm=18)*3; + + + + diff --git a/tel/TPmatlab/initiation/somme1.m b/tel/TPmatlab/initiation/somme1.m new file mode 100644 index 0000000..a68fa27 --- /dev/null +++ b/tel/TPmatlab/initiation/somme1.m @@ -0,0 +1,10 @@ +function x=somme1(n,m) +% déterminer ce que fait cette fonction à deux variables d'entrée. + +for i= 1:n + xtemp=0; + for j = 1:m + xtemp=xtemp+log(i)*exp(-j^2); + end + x(i)=xtemp; +end diff --git a/tel/TPmatlab/integration/TP3a/int_fcn.m b/tel/TPmatlab/integration/TP3a/int_fcn.m new file mode 100644 index 0000000..bc8fd10 --- /dev/null +++ b/tel/TPmatlab/integration/TP3a/int_fcn.m @@ -0,0 +1,72 @@ +function res=int_fcn(c,N,A,B,fcn,varargin) + + +% int_fcn : calcul approché d'intégrale par méthode classique +% +% ********************************************************* +% +% res=int_fcn(c,N,A,B,fcn,varargin) renvoie la valeur approchée de l'intégrale +% d'une fonction sur un intevalle par la méthodes des rectangles, des milieux, +% des trapèzes ou de simpson. +% +% variables d'entrées : +% * c définit la méthode : +% c=1 : méthode des rectangles, +% c=2 : méthode des trapèzes, +% c=3 : méthode des milieux, +% c=4 : méthode de Simpson; +% * N est le nombre de sous intervalles utilisés (correspondant à h=(B-A)/N) +% * [A,B] est l'intervalle d'intégration +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% * varargin sont les arguments optionnels de la fonction. +% attention, intégration par rapport à la première variable de f. +% +% variables de sortie +% * res est la valeur approchée de l'intégrale calculée. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + +% Contrôles d'entrée +% nombre d'arguments +if (A>=B) + error('il faut A=1) + P(2,1:2)=[1 0]; + Lb=1; + La=[1 0]; + for k=1:n-1 + Lc=[2*La 0]-[0 0 Lb]; + Lb=La; + La=Lc; + P(k+2,1:k+2)=Lc; + end + end +end diff --git a/tel/TPmatlab/integration/TP3b/demo_trace_poly_gauss.m b/tel/TPmatlab/integration/TP3b/demo_trace_poly_gauss.m new file mode 100644 index 0000000..1404755 --- /dev/null +++ b/tel/TPmatlab/integration/TP3b/demo_trace_poly_gauss.m @@ -0,0 +1,64 @@ +% script demo_trace_poly_gauss : tracé des polynômes de gauss P_0, P_1,...,P_n +% +% +% pour l'une des quatre familes : +% Gauss-Legendre, +% Gauss-Tchebytchev, +% Gauss-Hermite, +% Gauss-Laguerre. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% calcul_poly_gauss +% +% ********************************************************* + + +n=input('entrez l''entier naturel n : '); +N=input('entrez le nombre de points pour le graphique : '); +disp('choix de la famille : '); +disp('c=1 : Gauss-Legendre'); +disp('c=2 : Gauss-Tchebytchev'); +disp('c=3 : Gauss-Hermite'); +disp('c=4 : Gauss-Laguerre'); +c=input('entrez c : '); +if (c==3) | (c==4) + A=input('entrez la borne supérieure de l''intervalle d''etude : '); +else + A=1; +end +choix=input('voulez vous voir les légendes ? entrez o ou n : ','s'); +if (c==4) + h=A/(N); + x=0:h:A; +else + h=2*A/(N); + x=-A:h:A; +end +P=calcul_poly_gauss(c,n); +clf; +hold on; +if strcmp(choix,'o') + touslesnoms=[]; + toutesleslegendes=[]; + Q=5; + Qtestt=0; + for k=0:n + nom =num2str(k); + long=length(nom); + Qtestt=Qtestt+1; + touslesnoms=[touslesnoms,' ',nom]; + legende=[nom,blanks(Q-long+1)]; + toutesleslegendes=[toutesleslegendes;legende]; + end +end +for k=0:n + cou=couleur7(k); + plot(x,polyval(P(k+1,1:k+1),x),cou); +end +if strcmp(choix,'o') + legend(strvcat(toutesleslegendes),0); +end +hold off; diff --git a/tel/TPmatlab/integration/TP3b/int_gauss.m b/tel/TPmatlab/integration/TP3b/int_gauss.m new file mode 100644 index 0000000..63682f8 --- /dev/null +++ b/tel/TPmatlab/integration/TP3b/int_gauss.m @@ -0,0 +1,60 @@ +function Iapp=int_gauss(n,points,poids,rcn) + +% int_gauss : calcul de l'intégrale par méthode de Gauss à n+1 points +% +% ********************************************************* +% +% Iapp=int_gauss(n,points,poids,fcn) : valeur approchée de l'intégrale +% de la fonction rcn=r=f/w (w poids poids) sur l'intervalle sous la forme conventionnelle (a,b) avec +% l'une des quatre méthodes à n+1 points: +% Legendre, Tchebytchev, Hermite ou Laguerre, éventuellement en formel +% +% +% variables d'entrées : +% * n : entier naturel. +% * points : tableau de taille (q+1)*(q+1) (avec n<=q) tel que +% pour tout l dans {0,...,q}, points(l+1,1:l+1) contient +% les l+1 points x_i de la formule (à l+1 points) +% * poids : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient +% les l+1 poids D_i de la formule (à l+1 points) +% * rcn est une chaîne de caractère représentant la fonction r=f/w +% (de type inline, builtin ou par fichier M-file); +% variables de sortie +% * Iapp : valeur approchée de l'intégrale de f sur un intervalle conventionnel. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=4 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end +[q1,q2]=size(points); +[q3,q4]=size(poids); +if (q1~=q2) + error('attention, le troisième argument doit être une matrice carrée'); +end +if (q3~=q4) + error('attention, le quatrième argument doit être une matrice carrée'); +end +if (q1~=q3) + error('attention, le troisième et le quatrième argument doivent avoir la même taille'); +end +if (n>q1-1) + error('attention, n est trop grand'); +end + +% Corps d'algorithme. +Iapp=sum(feval(rcn,(points(n+1,1:n+1))).*poids(n+1,1:n+1)); diff --git a/tel/TPmatlab/integration/TP3b/int_gauss_legendre.m b/tel/TPmatlab/integration/TP3b/int_gauss_legendre.m new file mode 100644 index 0000000..dae77c4 --- /dev/null +++ b/tel/TPmatlab/integration/TP3b/int_gauss_legendre.m @@ -0,0 +1,63 @@ +function Iapp=int_gauss_legendre(n,a,b,points,poids,fcn) + +% int_gauss_legendre : calcul de l'intégrale par méthode de Gauss-Legendre à n+1 points +% +% ********************************************************* +% +% Iapp=int_gauss_legendre(n,a,b,points,poids,fcn) : valeur approchée de l'intégrale +% de la fonction fcn sur l'intervalle [a,b] avec la méthode de Gauss-Legendre +% à n+1 points. +% +% variables d'entrées : +% * n : entier naturel. +% * a,b : bornes de l'intervalle +% * points : tableau de taille (q+1)*(q+1) (avec n<=q) tel que +% pour tout l dans {0,...,q}, points(l+1,1:l+1) contient +% les l+1 points x_i de la formule (à l+1 points) +% * poids : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient +% les l+1 poids D_i de la formule (à l+1 points) +% * fcn est une chaîne de caractère représentant la fonction +% (de type inline, builtin ou par fichier M-file); +% variables de sortie +% * Iapp : valeur approchée de l'intégrale +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% +% + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=6 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end +if (a>=b) + error('attention, a doit être plus petit que b'); +end +[q1,q2]=size(points); +[q3,q4]=size(poids); +if (q1~=q2) + error('attention, le troisième argument doit être une matrice carrée'); +end +if (q3~=q4) + error('attention, le quatrième argument doit être une matrice carrée'); +end +if (q1~=q3) + error('attention, le troisième et le quatrième argument doivent avoir la même taille'); +end +if (n>q1-1) + error('attention, n est trop grand'); +end + +% Corps d'algorithme. +auxi=(b-a)/2; +Iapp=auxi*sum(feval(fcn,auxi*(points(n+1,1:n+1)+1)+a).*poids(n+1,1:n+1)); diff --git a/tel/TPmatlab/integration/TP3b/points_poids_gauss.m b/tel/TPmatlab/integration/TP3b/points_poids_gauss.m new file mode 100644 index 0000000..4ab79b3 --- /dev/null +++ b/tel/TPmatlab/integration/TP3b/points_poids_gauss.m @@ -0,0 +1,72 @@ +function [points,poids]=points_poids_gauss(c,n) + + +% points_poids_gauss : calcul des poids et des points de Gauss +% +% ********************************************************* +% +% [points,poids]=points_poids_gauss(c,n) renvoie les points D_i et les points x_i +% des formules d'intégration à l+1 points pour 0<=l<=n, pour +% l'une des quatre méthodes : Legendre, Tchebytchev, Hermite ou Laguerre. +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : méthode de Gauss-Legendre, +% c=2 : méthode de Gauss-Tchebychev, +% c=3 : méthode de Gauss-Hermite, +% c=4 : méthode de Gauss-Laguerre; +% +% variables de sortie +% * points : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, points(l+1,1:l+1) contient +% les l+1 points x_i de la formule (à l+1 points) +% * poids : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient +% les l+1 poids D_i de la formule (à l+1 points) +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% calcul_poly_gauss +% +% ********************************************************* + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end +if ~((c==1) | (c==2) | (c==3) | (c==4)) + error('c doit être un entier égal à 1, 2, 3 ou 4'); +end + +% Corps d'algorithme. +poids=zeros(n+1,n+1); +points=zeros(n+1,n+1); +if (c~=2) + P=calcul_poly_gauss(c,n+1); + for k=0:n + C=(sort(roots(P(k+2,1:k+2))))'; + [points(k+1,1:k+1),ind]=sort(C); + switch c + case 1 + auxipoids=2./((k+1)*polyval(polyder(P(k+2,1:k+2)),C).*polyval(P(k+1,1:k+1),C)); + case 3 + auxipoids=(2^(k+2)*sqrt(pi)*factorial(k+1))./((polyval(polyder(P(k+2,1:k+2)),C)).^2); + case 4 + auxipoids=1./(C.*(polyval(polyder(P(k+2,1:k+2)),C)).^2); + end + poids(k+1,1:k+1)=auxipoids(ind); + end +else + for k=0:n + points(k+1,1:k+1)=cos(((2*(0:k)+1)/(k+1))*pi/2); + poids(k+1,1:k+1)=(pi/(k+1))*ones(1,k+1); + end +end diff --git a/tel/TPmatlab/integration/TP3b/pointspoids_manuel.mat b/tel/TPmatlab/integration/TP3b/pointspoids_manuel.mat new file mode 100644 index 0000000..f2001c0 Binary files /dev/null and b/tel/TPmatlab/integration/TP3b/pointspoids_manuel.mat differ diff --git a/tel/TPmatlab/integration/TP3c/demo_tp3C.m b/tel/TPmatlab/integration/TP3c/demo_tp3C.m new file mode 100644 index 0000000..c73e9e8 --- /dev/null +++ b/tel/TPmatlab/integration/TP3c/demo_tp3C.m @@ -0,0 +1,36 @@ +% fichier de démo du tp3C +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +% gener_signal,determin_periode,visu_integ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +global etendue +etendue=1.2; +% valeur modifiable mais attention aux effets éventuels; la méthode +% est pensée pour cette valeur, conforme au problème initial. + + +% génération d'un signal pseudo-aléatoire sinusoidal +% ou acquisition de véritables mesures mises au format voulu +[X,Z,T_th]=gener_signal(4,1000); + +% recherche de période +[per,ind]=determin_periode(X,Z); + +% préparation de visualisation pour un 17 pouces +% vous pouvez adapter en créant une figure, en récupérant son champ +% position tel que vous le souhaitez; passez alors ces données +% dans l'affectation de fen_pos ci-dessous. +% close all; +% fen_pos=[385 45 620 640]; +% figure(1);set(1,'position',fen_pos); + +% Vous pouvez aussi lancer l'exécution une fois, mettre la fenetre +% à la taille souhaitée et avant une nouvelle exécution taper clf. + +% visualisation animée +visu_integ(X,Z,ind); + +% variantes utilisateur! + diff --git a/tel/TPmatlab/integration/TP3c/determin_periode.m b/tel/TPmatlab/integration/TP3c/determin_periode.m new file mode 100644 index 0000000..bc2d99b --- /dev/null +++ b/tel/TPmatlab/integration/TP3c/determin_periode.m @@ -0,0 +1,90 @@ +function [per,num]=determin_periode(X,Z) +% détermine la période d'un signal sinusoidal par calcul d'intégrale +% +% variables d'entrée +% X désigne le vecteur des dates régulièrement réparties du signal +% Z désigne le vecteur des images du signal s : Z=s(X). +% Il est essentiel que le signal issu de mesures a été toiletté +% préalablement sinon les calculs d'intégrales menés par exemple +% n'ont pas de sens. +% +% variable de sortie +% per est la valeur approchée de la période obtenue +% grace à la méthode retenue: voir sujet et fichier .doc. +% num désigne l'indice du vecteur X tel que X(num)-X(1)=per. +% On pourrait renvoyer comme paramètre de sortie aussi les deux +% premiers indices en lesquels Z est nul. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% simpson, +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Tests d'entrée +if nargin~=2 + error('Nombre de champs incorrect'); +end +% controle des données +if size(X)~=size(Z) + error('Champs incorrects. Voir controle des données'); +end + + +format long e; +% + + +% paramètres indispensables +pas_X=X(2)-X(1); +lgr=size(X,2); + +% détermination du deuxième changement de signe de Z +positifs=(Z>=0); +chgt=find(diff(positifs)); + +% Initialisation du cacul de l'intégrale de z=f(x) +% entre X0 et le deuxième changement de signe de Z. +% ind_deb, choisi pour disposer d'un nombre de points impairs, +% sera l'indice de début dans le calcul des intégrales suivantes. +ind_deb=chgt(2)+mod(chgt(2)+1,2); +integ0=simpson(pas_X,Z,1,ind_deb); +ind_fin=lgr-mod(lgr-ind_deb,2); +integ=integ0; + +% recherche dichotomique de l'intégrale nulle +% en fin de boucle +% l'integrale entre 1 et ind_deb est d'un signe +% l'integrale entre 1 et ind_fin est de l'autre +% donc l'intégrale entre 1 et leur milieu est considérée nulle. + +while (ind_fin-ind_deb>2) + mil=(ind_deb+ind_fin)/2; + % vérification de la conformité des indices deb et fin + % assurant un nombre impair de points pour intégrer. + mil=mil-mod(mil-ind_deb,2); + integ_ajout=simpson(pas_X,Z,ind_deb,mil); + if (integ0*(integ+integ_ajout)>=0) + ind_deb=mil; + integ=integ+integ_ajout; + else + ind_fin=mil; + end +end + +% sortie des paramètres attendus +num=(ind_deb+ind_fin)/2; +per=(num-1)*pas_X; +% on pourrait sortit un vecteur ind_zeros_Z qui contiendrait +% [chgt(1) chgt(2)]. + + + + + + + + + + + + + diff --git a/tel/TPmatlab/integration/TP3c/gener_signal.m b/tel/TPmatlab/integration/TP3c/gener_signal.m new file mode 100644 index 0000000..90288ba --- /dev/null +++ b/tel/TPmatlab/integration/TP3c/gener_signal.m @@ -0,0 +1,75 @@ +function [X,Z,T_th]=gener_signal(ampl,nb_int) +% génère un signal sinusoidal aléatoire +% +% variables d'entrée +% ampl désigne l'amplitude maximale du signal; +% nb_int désigne le nombre d'intervalles entre mesures; +% ainsi le vecteur des abscisses X comporte nb_int+1 valeurs. +% Dans cette application on le minore à 500, en réalité +% on acquiert 10000 mesures en 0.8 secondes environ! +% +% variables de sortie +% X est le vecteur des abscisses +% réparties sur 1,2 période environ comme dans l'application +% industrielle qui a fondé ce développement; ce paramètre sera +% global aux blocs qui le nécessitent sous le nom etendue. +% Z est le vecteur des images associé à X défini par +% Z= a*sin(omega*x+phi) +% où a,omega et phi sont pseudo-aléatoires, définis plus loin. +% T_th désigne la valeur théorique de période, sortie à des fins de +% vérification. Elle n'a pas de sens dans l'application industrielle +% puisqu'elle n'est pas connue! +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Tests d'entrée +if nargin~=2 + error('Nombre de champs incorrects'); +end +if (ampl>0)*(nb_int>500)==0 + error('Champs incorrects. Voir tests d''entrée'); +end + +global etendue +% valeur modifiable mais attention aux effets éventuels; la méthode +% est pensée pour cette valeur, conforme au problème initial. + +format long e; + +% Génération du signal "aléatoire" +% remplacé dans le réel par les données des capteurs +% renvoyées au format matlab choisi ci-dessous. +% Si t est élement de [0,1] comme les renvoie rand, alors +% T=(b-a)t+a est élément de [a,b]. + +aleas=rand(1,3); +if aleas(2)==0 + aleas(2)=1; +end + +% définition de l'amplitude a choisie entre ampl/2 et ampl, +% pour assurer la visibilité du signal produit. +a=ampl/2*aleas(1)+ampl/2; +% détermination de phi et omega +omega=abs(aleas(2)); +% période théorique +T_th=2*pi/omega; +% phi est ramené dans [O,T_th] +phi=T_th*aleas(3); + +% fabrication du vecteur X +pas_X=etendue*T_th/nb_int; +i=0:nb_int; +X=pas_X.*i; + +% fabrication de Z +Z=a.*sin(omega.*X+phi); + + + + + + diff --git a/tel/TPmatlab/integration/TP3c/simpson.m b/tel/TPmatlab/integration/TP3c/simpson.m new file mode 100644 index 0000000..2b39aec --- /dev/null +++ b/tel/TPmatlab/integration/TP3c/simpson.m @@ -0,0 +1,40 @@ +function integ=simpson(h,y,ind_deb,ind_fin) +% calcule de façon simplifiée une intégrale en méthode de simpson +% +% variables d'entrée +% h est le demi-pas d'intégration au sens du cours; +% y est le vecteur des valeurs de la fonction à intégrer; +% ind_deb est l'indice de la borne basse de l'intervalle d'intégration; +% ind_fin est l'indice de la borne haute de l'intervalle d'intégration; +% Il est implicite que le nombre (ind_fin-ind_deb) est un entier pair +% c'est-à-dire que le nombre de points d'intégration est impair. +% +% variables de sortie +% integ est la valeur approchée par Simpson de l'intégrale cherchée +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% Tests d'entrée +if nargin~=4 + error('Nombre de champs incorrect'); +end +% controle de données +bool1=(h>0)*(mod(ind_fin-ind_deb,2)==0)*(ind_deb0 + error('Champs incorrects. Voir controle des champs.'); +end + +global etendue +% Voir remarques dans gener_signal ou demo_tp3C + +% paramètres généraux +lgr=size(X,2); +pas_X=X(2)-X(1); + +% calibrages + +% pour les fenetres 1 et 2. +ampl1=etendue*max(Z); +lim1=[X(1) X(lgr) -ampl1 ampl1]; + +% outil intermédiaire : calcul de l'intégrale sur une demi-période +% détermination du deuxième changement de signe de Z. +% Peut etre remplacé par un passage de ces données issues +% de determin_periode, ce qui éviterait leur calcul ici. +positifs=(Z>=0); +chgt=find(diff(positifs)); + +ampl3=etendue*... + abs(simpson(pas_X,Z,chgt(1),chgt(2)+mod(chgt(2)-chgt(1),2))); +lim3=[-ampl3 ampl3 -3 4]; + + +% tracé du signal brut +subplot(3,1,1); +plot(X,Z); +axis(lim1); +legend('Signal acquis',0); + +% tracé de l'intégrale en cours de calcul + +% choix de la finesse de découpe de l'intégrale à calculer +% à modifier éventuellement. +nb_iterations=50; unite=10/nb_iterations; +pas=floor((lgr-chgt(2))/nb_iterations); + +% initialisation pour la deuxième fenetre +subplot(3,1,2); +title('Evolution de l''intégrale du signal'); +deb=1;fin=chgt(2)+mod(chgt(2)-deb,2); +area(X(deb:fin),Z(deb:fin)); +axis(lim1); +hold on; + +% initialisation de la troisieme fenetre +h_3=subplot(3,1,3); +set(h_3,'YTick',[]); +title('Evolution de la valeur de l''intégrale'); +integ0=simpson(pas_X,Z,deb,fin); +area([0 integ0],[1 1]);axis(lim3); +pause(0.5); +hold on; + +% initialisations avant boucle +deb=fin;fin=deb+pas+mod(pas,2); +ajout=simpson(pas_X,Z,deb,fin); +integ=integ0+ajout; + +while integ*integ0>=0 + + % visualisation de l'intégrale en cours de calcul + rx2=X(deb:fin);rz2=Z(deb:fin); + subplot(3,1,2); + title('Evolution de l''intégrale du signal'); + area(rx2,rz2);axis(lim1); + hold on; + + % visualisation de la valeur de l'intégrale + h_3=subplot(3,1,3); + set(h_3,'YTick',[]); + title('Evolution de la valeur de l''intégrale'); + cla; + area([0 integ],[1 1]);axis(lim3); + + % préparation de boucle suivante + pause(unite); + deb=fin;fin=deb+pas+mod(pas,2); + ajout=simpson(pas_X,Z,deb,fin); + integ=integ+ajout; + + % traitement particulier de fin de boucle + if fin>ind_int_nulle + + % fenetre 2 finale + rx2=X(1:ind_int_nulle);rz2=Z(1:ind_int_nulle); + subplot(3,1,2); + title('Evolution de l''intégrale du signal'); + cla; + area(rx2,rz2);axis(lim1); + % fenetre 3 finale + h_3=subplot(3,1,3); + set(h_3,'YTick',[]); + title('Evolution de la valeur de l''intégrale'); + cla; + area([0 ampl3/100],[1 1]);axis(lim3); + end + +end diff --git a/tel/TPmatlab/integration/TP3d/choix_du_support.m b/tel/TPmatlab/integration/TP3d/choix_du_support.m new file mode 100644 index 0000000..1b0ab07 --- /dev/null +++ b/tel/TPmatlab/integration/TP3d/choix_du_support.m @@ -0,0 +1,72 @@ +function [W1,W2,W3]=choix_du_support(pd1,pd2,pd3,v1,v2,v3,verif) +% calcule les poids explicites d'une formule de quadrature et vérifie +% son caractère exact pour les fonctions polynomes de P2. +% +% variables d'entrée +% pd1,pd2,pd3 sont les symboliques issus de la fonction +% recherche_poids pour le problème considéré. +% v1,v2,v3 sont trois valeurs distinctes de l'intervalle [-1,1] +% qui particularisent les valeurs formelles. +% verif est un booléen +% si verif=1: le logiciel vérifie formellement que les intégrales +% exacte et approchée coincident pour les fonctions considérées; +% sinon, il omet cette vérification. +% +% variables de sortie +% W1,W2,W3 désignent les poids de la formule de quadrature. +% Ainsi l'intégrale I de la fonction f mesurée est donnée par: +% I=W1f(v1)+W2f(v2)+W3f(v3). +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% pour l'exemple, intervention de globales! +global t1 t2 t3; + +% tests d'entrée +if nargin~=7 + error('Nombre des arguments incorrect'); +end +v=diff(sort([v1 v2 v3])); +if max(v.^2)>1 + error('Les points ne sont pas tous dans [-1,1].Voir controle.'); +end +if sum(v==0)>0 + error('Les points ne sont pas distincts. Voir controle.'); +end + +% algorithme proprement dit +v=[v1 v2 v3]; +v=sym(v,'r'); +poids=[pd1 pd2 pd3]; +poids_expl=sym(subs(poids,[t1 t2 t3],[v1 v2 v3]),'r'); +W1=poids_expl(1);W2=poids_expl(2);W3=poids_expl(3); + +% traitement de la vérification si l'utilisateur l'a demandée + +if verif==1 + syms a b c x; + ch1='On veut intégrer sur [-1,1] la fonction f donnée par'; + ch2='f(x)=ax^2+bx+c, où a,b,c sont symboliques.'; + ch3='La valeur exacte de l''intégrale est :'; + ch4='La valeur fournie par formule de quadrature est :'; + ch5='obtenue à partir de la valeur du signal aux dates'; + ch6='avec les poids respectifs'; + ch7='L''écart entre les deux expressions est de :'; + fx=a*x^2+b*x+c; + v_app=W1*subs(fx,x,v1)+W2*subs(fx,x,v2)+W3*subs(fx,x,v3); + % qu'on peut remplacer par un produit matriciel! + v_ex=int(fx,x,-1,1); + % l'indication de x est inutile en fait... + ecart=v_app-v_ex; + + % affichage des résultats + disp(ch1);disp(ch2); + disp(ch3);v_ex + disp(ch4);v_app + disp(ch5);v + disp(ch6);poids_expl + disp(ch7);ecart +end + diff --git a/tel/TPmatlab/integration/TP3d/demo_tp3D.m b/tel/TPmatlab/integration/TP3d/demo_tp3D.m new file mode 100644 index 0000000..d05cae2 --- /dev/null +++ b/tel/TPmatlab/integration/TP3d/demo_tp3D.m @@ -0,0 +1,26 @@ +% fichier de démo du tp3D +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +% recherche_poids,choix_du_support +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +global t1 t2 t3; +syms t1 t2 t3; + +% solution du problème formel +[w1,w2,w3]=recherche_poids(t1,t2,t3); + +% particularisation des valeurs des dates retenues + +% pour demande de vérification formelle +verif=1; + +% variante 1 +[W1,W2,W3]=choix_du_support(w1,w2,w3,-0.5,0,0.5,verif); + +% variante 2 +% [W1,W2,W3]=choix_du_support(w1,w2,w3,-sqrt(2)/2,0,sqrt(2)/2,verif); + +% variantes utilisateur! + diff --git a/tel/TPmatlab/integration/TP3d/recherche_poids.m b/tel/TPmatlab/integration/TP3d/recherche_poids.m new file mode 100644 index 0000000..9cd8f53 --- /dev/null +++ b/tel/TPmatlab/integration/TP3d/recherche_poids.m @@ -0,0 +1,31 @@ +function [w1,w2,w3]=recherche_poids(t1,t2,t3) +% recherche les poids d'une formule de quadrature sur [-1,1] +% exacte pour toute les fonctions polynomes de degré au plus deux; +% Le traitement est effectué en symbolique. +% +% variables d'entrée +% t1,t2,t3 désignent trois dates distinctes de [-1,1] en lesquelles +% sont effectuées les mesures de la fonction à intégrer; +% elles serviront de support pour la formule de quadrature. +% +% variables de sortie +% w1,w2,w3 sont les poids associés qui rendent l'intégration exacte +% pour toute fonction polynome de degré au plus deux. Voir documentation +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% tests d'entrée +if nargin~=3 + error('Nombre d''arguments incorrect.'); +end + +% corps d'algorithme proprement dit +modele=[t1,t1^2]'; +mat=sym(ones(3)); +mat(2:3,1)=modele; +mat(2:3,2)=subs(modele,t1,t2);mat(2:3,3)=subs(modele,t1,t3); +sol=(mat\sym([2 0 2/3]'))'; +w1=sol(1);w2=sol(2);w3=sol(3); + diff --git a/tel/TPmatlab/integration/TP3e/base_pn.m b/tel/TPmatlab/integration/TP3e/base_pn.m new file mode 100644 index 0000000..359639a --- /dev/null +++ b/tel/TPmatlab/integration/TP3e/base_pn.m @@ -0,0 +1,36 @@ +function[base]=base_pn(x,X) +% engendre la base de Pn utilisée dans l'écriture de Newton +% à partir des symboliques contenus dans x. +% +% variables d'entrée +% x est un vecteur de deg_max+1 symboliques repréentant les points +% du support d'interpolation. +% X représente la variable symbolique +% générique dans les écritures polynomiales. +% +% variables de sortie +% base est un vecteur de symboliques représentant la base +% 1,(X-x0),...,(X-x0)...(X-xn-1) de Pn. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% tests d'entrée +if nargin~=2 + error('Nombre d''arguments incorrect.'); +end + +% corps d'algorithme proprement dit + +% évite de passer une globale +deg_max=size(x,2)-1; + +base=sym(ones(1,deg_max+1));gen=X*ones(1,deg_max)-x(1:deg_max); +for i=2:deg_max+1, + base(i)=prod(gen(1:i-1)); +end + + + diff --git a/tel/TPmatlab/integration/TP3e/demo_tp3E.m b/tel/TPmatlab/integration/TP3e/demo_tp3E.m new file mode 100644 index 0000000..d0717c0 --- /dev/null +++ b/tel/TPmatlab/integration/TP3e/demo_tp3E.m @@ -0,0 +1,66 @@ +% fichier script de demo_tp3E +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +% outils_pol (donc diff_div_dist et base_pn),deriv_app +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% initialisations diverses +global d base; + +% indéterminée classique des polynomes +syms X real; +deg_max=6; + +% Création des variables symboliques nécessaires; +% pourrait se faire dans une fonction à part! +x=[];y=[]; +ch1='x';ch2='y'; +for i=0:deg_max + sym([ch1 num2str(i)],'real'); + x=[x sym([ch1 num2str(i)])]; + sym([ch2 num2str(i)],'real'); + y=[y sym([ch2 num2str(i)])]; +end + +% détermination une seule fois des outils pour le calcul +% du polynome d'interpolation maximal. +[d,base]=outils_pol(x,y,X); + +% boucle des recherches de formules de dérivation approchée + +% initialisation +rep=1; +while rep==1 + + % Choix de l'ordre de dérivation et du degré du polynome d'interpolation; + % les acharnés échapperont au controle! + k=-1;compt=0; + while ((k<=0)|(k>deg_max))&(compt<4) + ch1=input('Choisissez un ordre de dérivation k (0deg_max))&(compt<4) + ch2=input('Choisissez un degré n (n<7)de polynome d''interpolation cohérent\n','s'); + n=floor(str2num(ch2)); + compt=compt+1; + % si compt==4 message d'injure... + end + + % Calcul proprement dit + res=deriv_app(k,n,x,y,X); + + % Désir d'un autre calcul + ch=input('Voulez-vous effectuer d''autres calculs ?(si oui: taper 1) \n','s'); + rep=str2num(ch); +end + + + + + diff --git a/tel/TPmatlab/integration/TP3e/deriv_app.m b/tel/TPmatlab/integration/TP3e/deriv_app.m new file mode 100644 index 0000000..d39f93e --- /dev/null +++ b/tel/TPmatlab/integration/TP3e/deriv_app.m @@ -0,0 +1,81 @@ +function res=deriv_app(k,n,x,y,X) +% fournit une valeur approchée de dérivée d'ordre k avec un polynome de degré n +% +% variables d'entrée +% k est un entier naturel qui désigne l'ordre de dérivation demandé +% n est un entier désignant le degré du polynome d'interpolation utilisé pour +% la dérivation approchée. +% Selon la parité de n, +% f est supposée donnée en x=(x0,...,xn) centrés autour de t, +% par les valeurs y=(y0,...,yn) via yi=f(xi). +% X désigne l'indéterminée classique. +% +% variable de sortie +% le résultat res uniquement défini comme un symbolique. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +global d base; +syms t h real; + +% tests d'entrée +if nargin~=5 + error('Nombre des arguments incorrect'); +end + +% effectué en amont... On peut choisir un autre style. + +% algorithme proprement dit + + +% calcul du polynome interpolant adapté aux choix effectués +% et calcul de sa dérivée formelle +util=d(1:n+1)*(base(1:n+1))'; +deriv=diff(util,'X',k); + +% détermination du support centré selon parité de n +% sachant qu'il y a n+1 points au support. +if (2*fix(n/2)==n) + % n est pair + % détermination du centre parmi les valeurs xi + c=n/2+1; + + % écriture des centres + xcentre(c)=t; + for i=n/2:-1:1, + xcentre(c-i)=t-i*h;xcentre(c+i)=t+i*h; + end +else + % n est impair + sup=(n+1)/2; + for i=1:sup, + xcentre(i)=t-(sup-i+1)*h;xcentre(n+2-i)=t+(sup-i+1)*h; + end +end + +% valeur approchée de la dérivée en fonction des valeurs centrées +res1=simple(subs(deriv,x(1:n+1),xcentre)); + +% sortie des résultats +disp('Les valeurs considérées en x sont :'); +disp(x(1:n+1)); +disp('soit encore, une fois centrées :'); +disp(xcentre); +disp('auxquelles correspondent respectivement en f(x):'); +disp(y(1:n+1)); +sor=[]; +if (2*fix(n/2)==n) + sor=sprintf([sor '=D%df(t=x%d)'],k,c-1); +end +ch=sprintf(['La valeur approchée de D%df(t)' sor ' est,'],k); +disp(ch); +disp('sous forme classique :'); +pretty(simple(subs(res1,X,t))); +res=simple(subs(res1,X,t)); + + + + + diff --git a/tel/TPmatlab/integration/TP3e/outils_pol.m b/tel/TPmatlab/integration/TP3e/outils_pol.m new file mode 100644 index 0000000..2985251 --- /dev/null +++ b/tel/TPmatlab/integration/TP3e/outils_pol.m @@ -0,0 +1,35 @@ +function[d,base]=outils_pol(x,y,X) +% construit les outils de détermination du polynome d'interpolation +% d'une fonction donnée. +% +% variables d'entrée +% x est le vecteur symbolique des points de support d'interpolation; +% y est le vecteur symbolique des valeurs prises par la fonction en +% ces points. +% X est une variable symbolique,nom de l'indéterminée. +% +% variables de sortie +% d est le vecteur symbolique des différences divisées associées; +% base représente la base de Newton associée en symbolique. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +% +% diff_div_dist, base_pn +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% tests d'entrée +if nargin~=3 + error('Nombre d''arguments incorrect.'); +end +% tests éventuels de controle sur les tailles respectives +% des symboliques x et y. + +% corps d'algorithme proprement dit + +d=diff_div_dist(x,y); +base=base_pn(x,X); + + + + diff --git a/tel/TPmatlab/integration/TP3f/demo_compar_methode_int.m b/tel/TPmatlab/integration/TP3f/demo_compar_methode_int.m new file mode 100644 index 0000000..7b53449 --- /dev/null +++ b/tel/TPmatlab/integration/TP3f/demo_compar_methode_int.m @@ -0,0 +1,68 @@ +% script : demo_compar_methode_int compare les quatre méthodes d'intégrations +% - des rectangles , +% - du point milieu , +% - des trapèzes , +% - de Simpson . +% en comparant l'intégrale calculée avec l'intégrale exacte. +% en terme d'erreur commise. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% saisiefonction, int_fcn +% +% ********************************************************* + + +% ATTENTION : si erreur nulle, log(0)/log(10)=log10d0; +log10d0=-20; + +a=input('Entrez l''extrémité inférieure de [a,b] : '); +b=input('Entrez l''extrémité supérieure de [a,b] : '); +nmin=input('Entrez la puissance de dix du nombre minimum de point de discrétisation : '); +nmax=input('Entrez la puissance de dix du nombre maximum de point de discrétisation : '); +p=input('Entrez le nombre total de point de calculs : '); +disp(' '); +disp(' '); +disp('entrée de la fonction à intégrer'); +ch1=saisiefonction; +disp(' '); +disp(' '); +disp('entrée de sa primitive'); +ch2=saisiefonction; +x=zeros(1,p); +erreur=zeros(4,p); +nombrecalcul=zeros(4,p); +% calcul des erreurs pour les quatre méthodes +integraleexacte=feval(ch2,b)-feval(ch2,a); +compteur=0; +for k=logspace(nmin,nmax,p) + compteur=compteur+1; + for methode=1:4 + kf=fix(k); + x(compteur)=kf; + intva=int_fcn(methode,kf,a,b,ch1); + err=abs(intva-integraleexacte); + if (err==0) + erreur(methode,compteur)=log10d0; + else + erreur(methode,compteur)=log10(err); + end + end +end +log10x=log10(x); +% tracé des erreurs commises par les différentes méthodes +plot(log10x,erreur(1,:),'b',... + log10x,erreur(2,:),'r',... + log10x,erreur(3,:),'g',... + log10x,erreur(4,:),'y'); +if exist(ch1)==1 + title(['erreur d''intégration pour la fonction ','inline(',formula(ch1),')']); +else + title(['erreur d''intégration pour la fonction ',ch1]); +end +xlabel('logarithme 10 du nombre de points d''intégration'); +ylabel('logarithme 10 de l''erreur d''intégration'); +legend('rectangle','trapèzes','milieux','Simpson',0); + + diff --git a/tel/TPmatlab/integration/TP3f/demo_compar_methode_int_inc.m b/tel/TPmatlab/integration/TP3f/demo_compar_methode_int_inc.m new file mode 100644 index 0000000..84c65ec --- /dev/null +++ b/tel/TPmatlab/integration/TP3f/demo_compar_methode_int_inc.m @@ -0,0 +1,68 @@ +% script : demo_compar_methode_int_inc compare les quatre méthodes d'intégrations +% - des rectangles , +% - du point milieu , +% - des trapèzes , +% - de Simpson . +% en comparant +% en comparant l'intégrale calculée avec N points +% avec l'intégrale calculée avec 2N points +% en terme d'erreur. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% saisiefonction, int_fcn +% +% ********************************************************* + + +% ATTENTION : si erreur nulle, log(0)/log(10)=log10d0; +log10d0=-20; + +a=input('Entrez l''extrémité inférieure de [a,b] : '); +b=input('Entrez l''extrémité supérieure de [a,b] : '); +nmin=input('Entrez la puissance de dix du nombre minimum de point de discrétisation : '); +nmax=input('Entrez la puissance de dix du nombre maximum de point de discrétisation : '); +p=input('Entrez le nombre total de point de calculs : '); +disp(' '); +disp('entrée de la fonction à intégrer'); +ch1=saisiefonction; +x=zeros(1,p); +erreur=zeros(4,p); +temps=zeros(4,p); +nombrecalcul=zeros(4,p); +% calcul des erreurs pour les quatre méthodes +compteur=0; +for k=logspace(nmin,nmax,p) + compteur=compteur+1; + for methode=1:4 + kf=fix(k); + x(compteur)=kf; + intva=int_fcn(methode,kf,a,b,ch1); + intvab=int_fcn(methode,2*kf,a,b,ch1); + err=abs(intva-intvab); + if (err==0) + erreur(methode,compteur)=log10d0; + else + erreur(methode,compteur)=log10(err); + end + end +end + + +log10x=log10(x); +% tracé des erreurs commises par les différentes méthodes +plot(log10x,erreur(1,:),'b',... + log10x,erreur(2,:),'r',... + log10x,erreur(3,:),'g',... + log10x,erreur(4,:),'y'); +if exist(ch1)==1 + title(['erreur d''intégration pour la fonction ','inline(',formula(ch1),')']); +else + title(['erreur d''intégration pour la fonction ',ch1]); +end +xlabel('logarithme 10 du nombre de points d''intégration'); +ylabel('logarithme 10 de l''erreur d''intégration'); +legend('rectangle','trapèzes','milieux','Simpson',0); + + diff --git a/tel/TPmatlab/integration/TP3g/calcul_poly_gauss_symb.m b/tel/TPmatlab/integration/TP3g/calcul_poly_gauss_symb.m new file mode 100644 index 0000000..4c68424 --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/calcul_poly_gauss_symb.m @@ -0,0 +1,86 @@ +function P=calcul_poly_gauss_symb(c,n) + +% calcul_poly_gauss_symb : calcul des polynômes de Gauss en symbolique (récurrence) +% +% ********************************************************* +% +% P=calcul_poly_gauss_symb(c,n) renvoie les polynômes de Gauss P_0,P_1,...,P_n +% pour les quatre méthodes génériques (Legendre, Tchebytchev, Hermite +% et Laguerre; calcul par récurrence. +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : Gauss-Legendre, +% c=2 : Gauss-Tchebytchev, +% c=3 : Gauss-Hermite, +% c=4 : Gauss-Laguerre. +% variables de sortie +% * P : tableau de taille (n+1)*(n+1) tel que +% pour tout i dans {1,...,n+1}, L(i,1:i) contient +% le polynôme L_{i-1} (coefficient dominant en tête) +% à coefficients symboliques +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end +if ~((c==1) | (c==2) | (c==3) | (c==4)) + error('c doit être un entier égal à 1, 2, 3 ou 4'); +end + +% Corps d'algorithme. +P=zeros(n+1,n+1); +P=sym(P); +if (c~=2) + P(1,1)=1; + Lb=[]; + La=[1]; + for k=0:n-1 + sk=sym(k); + switch(c) + case 1 + auxi0=1/(sk+1); + auxi1=(2*sk+1)*auxi0; + auxi2=-sk*auxi0; + Lc=[auxi1*La 0]+[0 0 auxi2*Lb]; + case 3 + Lc=[2*La 0]-[0 0 2*sk*Lb]; + case 4 + auxi0=1/(sk+1); + auxi1=(2*sk+1)*auxi0; + auxi2=-sk*auxi0; + Lc=[-auxi0*La 0]+[0 auxi1*La]+[0 0 auxi2*Lb]; + end + Lb=La; + La=Lc; + P(k+2,1:k+2)=Lc; + end +else + P(1,1)=1; + if (n>=1) + P(2,1:2)=[1 0]; + Lb=1; + La=[1 0]; + for k=1:n-1 + Lc=[2*La 0]-[0 0 Lb]; + Lb=La; + La=Lc; + P(k+2,1:k+2)=Lc; + end + end +end diff --git a/tel/TPmatlab/integration/TP3g/calcul_poly_gauss_symb_dir.m b/tel/TPmatlab/integration/TP3g/calcul_poly_gauss_symb_dir.m new file mode 100644 index 0000000..cc93f5c --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/calcul_poly_gauss_symb_dir.m @@ -0,0 +1,85 @@ +function P=calcul_poly_gauss_symb_dir(c,n) + +% calcul_poly_gauss_symb_dir : calcul des polynômes de Gauss en symbolique directement +% +% ********************************************************* +% +% P=calcul_poly_gauss_symb_dir(c,n) renvoie les polynômes de Gauss P_0,P_1,...,P_n +% pour les quatre méthodes génériques (Legendre, Tchebytchev, Hermite +% et Laguerre; calcul par l'expression explicite. +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : Gauss-Legendre, +% c=2 : Gauss-Tchebytchev, +% c=3 : Gauss-Hermite, +% c=4 : Gauss-Laguerre. +% variables de sortie +% * P : tableau de taille (n+1)*(n+1) tel que +% pour tout i dans {1,...,n+1}, L(i,1:i) contient +% le polynôme L_{i-1} (coefficient dominant en tête) +% à coefficients symboliques +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end +if ~((c==1) | (c==2) | (c==3) | (c==4)) + error('c doit être un entier égal à 1, 2, 3 ou 4'); +end + +% Corps d'algorithme. +P=zeros(n+1,n+1); +P=sym(P); +syms x; +switch(c) + case 1 + P(1)=sym(1); + Q=sym(1); + a=sym(1); + b=sym(1); + for k=1:n + Q=Q*(x^2-1); + a=a*sym(2); + b=b*sym(k); + R=diff(Q,k); + P(k+1,1:k+1)=sym(sym2poly(expand(simplify(1/(a*b)*R)))); + end + case 2 + for k=0:n + P(k+1,1:k+1)=sym(sym2poly(expand(cos(sym(k)*acos(x))))); + end + case 3 + P(1)=sym(1); + R=exp(-x^2); + for k=1:n + R=diff(R); + P(k+1,1:k+1)=sym(sym2poly((-1)^k*expand(simplify(exp(x^2)*R)))); + end + case 4 + P(1)=sym(1); + a=sym(1); + c=sym(1); + for k=1:n + a=a*x; + c=c*sym(k); + Q=(exp(-x)*a); + R=diff(Q,k); + P(k+1,1:k+1)=sym(sym2poly(expand(simplify(exp(x)*R))/c)); + end +end + \ No newline at end of file diff --git a/tel/TPmatlab/integration/TP3g/demo_compar_methode_poly_gauss.m b/tel/TPmatlab/integration/TP3g/demo_compar_methode_poly_gauss.m new file mode 100644 index 0000000..21250ac --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/demo_compar_methode_poly_gauss.m @@ -0,0 +1,40 @@ +% demo_compar_methode_poly_gauss : ce script permet de comparer les deux méthodes +% de calcul des polynômes orthogonaux, +% (par récurrence et directement) en terme de temps de calcul en faisant un graphe +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% calcul_poly_gauss_symb, calcul_poly_gauss_symb_dir +% +% ********************************************************* + +clear all; +c=input('entrez c : '); +n=input('entrez n : '); +temps1=zeros(1,n+1); +nombreop1=zeros(1,n+1); +temps2=zeros(1,n+1); +nombreop2=zeros(1,n+1); +for k=0:n + tic; + d1=calcul_poly_gauss_symb_dir(c,k); + temps1(k+1)=toc; + tic; + d2=calcul_poly_gauss_symb(c,k); + temps2(k+1)=toc; +end +plot(1:n+1,temps1,1:n+1,temps2,'-.'); +legend('méthode directe','méthode par récurrence') +switch c + case 1 + title('polynômes de Gauss-Legendre'); + case 2 + title('polynômes de Gauss-Tchebychev'); + case 3 + title('polynômes de Gauss-Hermite'); + case 4 + title('polynômes de Gauss-Laguerre'); +end +xlabel('degré'); +ylabel('temps'); diff --git a/tel/TPmatlab/integration/TP3g/expression_part_symb.m b/tel/TPmatlab/integration/TP3g/expression_part_symb.m new file mode 100644 index 0000000..bc23eee --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/expression_part_symb.m @@ -0,0 +1,61 @@ +function res=expression_part_symb(c,n) + + +% expression_part_symb : calcul de l'intégrale de x^n*w(x) sur (a,b) (en symbolique) +% +% ********************************************************* +% +% expression_part_symb(c,n) : calcul en symbolique de \int_a^b x^n w(x) dx +% +% +% variables d'entrées : +% * n : entier naturel. +% * c définit le poids w : +% c=1 : Gauss-Legendre, +% c=2 : Gauss-Tchebytchev, +% c=3 : Gauss-Hermite, +% c=4 : Gauss-Laguerre. +% variables de sortie +% *res=valeur de l'inétrale +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if ~((c==1) | (c==2) | (c==3) | (c==4)) + error('c doit être un entier égal à 1, 2, 3 ou 4'); +end + +% Corps d'algorithme. +syms x a b; +switch c + case 1 + a=sym(-1); + b=sym(1); + w=sym(1); + case 2 + a=sym(-1); + b=sym(1); + w=1/sqrt(1-x^2); + case 3 + a=-sym(inf); + b=sym(inf); + w=exp(-x^2); + case 4 + a=sym(0); + b=sym(+inf); + w=exp(-x); +end +res=simplify(int(x^n*w,a,b)); diff --git a/tel/TPmatlab/integration/TP3g/points_poids_gauss_symb.m b/tel/TPmatlab/integration/TP3g/points_poids_gauss_symb.m new file mode 100644 index 0000000..98a5967 --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/points_poids_gauss_symb.m @@ -0,0 +1,106 @@ +function [points,poids]=points_poids_gauss_symb(c,n) + +% points_poids_gauss_symb : calcul des poids et des points de Gauss (symbolique) +% +% ********************************************************* +% +% [points,poids]=points_poids_gauss_symb(c,n) renvoie les points D_i et les points x_i +% des formules d'intégration à l+1 points pour 0<=l<=n, pour +% l'une des quatre méthodes : Legendre, Tchebytchev, Hermite ou Laguerre, +% en symbolique +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : méthode de Gauss-Legendre, +% c=2 : méthode de Gauss-Tchebychev, +% c=3 : méthode de Gauss-Hermite, +% c=4 : méthode de Gauss-Laguerre; +% +% variables de sortie +% * points : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, points(l+1,1:l+1) contient +% les l+1 points x_i de la formule (à l+1 points) +% * poids : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient +% les l+1 poids D_i de la formule (à l+1 points) +% +% Pour ces deux tableau, évaluation numérique (augmenter le nombre +% de chiffre avec la fonction digits) ou formelle si possible. +% méthode : points : recherche de zéro, +% poids : utilisation des formules du cours (dérivée) +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% calcul_poly_gauss_symb, polyval_symb, polyder_symb +% +% ********************************************************* + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end + +% Corps d'algorithme. +poids=zeros(n+1,n+1); +points=zeros(n+1,n+1); +poids=sym(poids); +points=sym(points); +spi=sym(pi); +if (c~=2) + L=calcul_poly_gauss_symb(c,n+1); + if ((1<=c) & (c<=3)) + points(1,1:1)=0; + else + points(1,1)=1; + end + switch c + case 1 + poids(1,1)=2; + case 3 + poids(1,1)=sqrt(spi); + case 4 + poids(1,1)=1; + end + if (c==3) + auxis=sym(4)*sqrt(spi); + end + syms x; + for k=1:n + sk=sym(k); + C=polyval_symb(L(k+2,1:k+2),x); + E=solve(C); + points(k+1,1:k+1)=E; + switch c + case 1 + Pp=polyder_symb(L(k+2,1:k+2)); + for j=1:k+1 + poids(k+1,j)=2./((sk+1)*polyval_symb(Pp,E(j)).*polyval_symb(L(k+1,1:k+1),E(j))); + end + case 3 + Pp=polyder_symb(L(k+2,1:k+2)); + auxis=2*auxis*(sk+1); + for j=1:k+1 + poids(k+1,j)=auxis/(polyval_symb(Pp,E(j)))^2; + end + case 4 + Pp=polyder_symb(L(k+2,1:k+2)); + for j=1:k+1 + poids(k+1,j)=1/(E(j)*(polyval_symb(Pp,E(j)))^2); + end + end + end +else + for k=0:n + sk=sym(k); + points(1+k,1:1+k)=cos(((2*sym(0:k))+1)/(sk+1)*spi/2); + poids(1+k,1:1+k)=spi/(sk+1)*ones(1,k+1); + end +end + diff --git a/tel/TPmatlab/integration/TP3g/points_poids_gauss_vand.m b/tel/TPmatlab/integration/TP3g/points_poids_gauss_vand.m new file mode 100644 index 0000000..e64b63b --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/points_poids_gauss_vand.m @@ -0,0 +1,95 @@ +function [points,poids]=points_poids_gauss_vand(c,n) + +% points_poids_gauss_vand : calcul des poids et des points de Gauss (avec Vandermonde) +% +% ********************************************************* +% +% [points,poids]=points_poids_gauss_vand(c,n) renvoie les points D_i et les points x_i +% des formules d'intégration à l+1 points pour 0<=l<=n, pour +% l'une des quatre méthodes : Legendre, Tchebytchev, Hermite ou Laguerre, +% éventuellement en formel +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : méthode de Gauss-Legendre, +% c=2 : méthode de Gauss-Tchebychev, +% c=3 : méthode de Gauss-Hermite, +% c=4 : méthode de Gauss-Laguerre; +% +% variables de sortie +% * points : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, points(l+1,1:l+1) contient +% les l+1 points x_i de la formule (à l+1 points) +% * poids : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient +% les l+1 poids D_i de la formule (à l+1 points) +% +% Pour ces deux tableau, évaluation numérique (augmenter le nombre +% de chiffre avec la fonction digits) ou formelle si possible. +% méthode : points : recherche de zéro, +% poids : calcul à partir des points en inversant la matrice de Vandermonde +% +% +% ************ Fonctions connexes utilisées *************** +% +% calcul_poly_gauss_symb, polyval_symb, expression_part_symb +% +% ********************************************************* + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end + +% Corps d'algorithme. +poids=zeros(n+1,n+1); +points=zeros(n+1,n+1); +poids=sym(poids); +points=sym(points); +spi=sym(pi); +if (c~=2) + L=calcul_poly_gauss_symb(c,n+1); + if ((1<=c) & (c<=3)) + points(1,1:1)=0; + else + points(1,1)=1; + end + switch c + case 1 + poids(1,1)=2; + case 3 + poids(1,1)=sqrt(spi); + case 4 + poids(1,1)=1; + end + U=zeros(n+1,1); + U=sym(U); + syms x; + for k=1:n + C=polyval_symb(L(k+2,1:k+2),x); + E=solve(C); + points(k+1,1:k+1)=E; + for i=0:k + U(i+1)=expression_part_symb(c,i); + end + M(1,1:k+1)=sym(ones(1,k+1)); + for i=2:1+k + for j=1:k+1 + M(i,j)=(points(k+1,j))^(i-1); + end + end + poids(k+1,1:k+1)=M\U(1:k+1); + end +else + for k=0:n + sk=sym(k); + points(1+k,1:1+k)=cos((2*sym(0:k)+1)/(sk+1)*spi/2); + poids(1+k,1:1+k)=spi/(sk+1)*ones(1,k+1); + end +end diff --git a/tel/TPmatlab/integration/TP3g/pointspoids_symb.mat b/tel/TPmatlab/integration/TP3g/pointspoids_symb.mat new file mode 100644 index 0000000..8ed31f4 Binary files /dev/null and b/tel/TPmatlab/integration/TP3g/pointspoids_symb.mat differ diff --git a/tel/TPmatlab/integration/TP3g/polyder_symb.m b/tel/TPmatlab/integration/TP3g/polyder_symb.m new file mode 100644 index 0000000..38b0d5e --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/polyder_symb.m @@ -0,0 +1,45 @@ +function res=polyder_symb(P) + +% polyder_symb : intégration formelle de polynôme +% +% +% ********************************************************* +% ATTENTION, polyder ne marche pas sur un tableau +% à coefficients formels !! +% +% ********************************************************* +% +% res=polyval_sym(P,x) : renvoie la dérivée de P en x (en formel) : +% si n+1 est la longueur de P : +% res=[n*P(1) (n-1)*P(2) ... 2*P(n-1) P(n)] +% +% variables d'entrées : +% * P : polynômes à coefficients formels +% variables de sortie : +% * res : tableau correspondant à la dérivée de P +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=1 + error('nombre d''arguments de la fonction incorrect'); +end + + +% Corps d'algorithme. +if (isempty(P)) | (P==0) + res=0; +else + n=length(P)-1; + res=[P(1:n).*sym(n:-1:1)]; +end \ No newline at end of file diff --git a/tel/TPmatlab/integration/TP3g/polyval_symb.m b/tel/TPmatlab/integration/TP3g/polyval_symb.m new file mode 100644 index 0000000..05e7e27 --- /dev/null +++ b/tel/TPmatlab/integration/TP3g/polyval_symb.m @@ -0,0 +1,48 @@ +function res=polyval_symb(P,x) + +% polyval_symb : évaluation formelle de polynôme +% +% ********************************************************* +% ATTENTION, ni polyval ni poly2sym ne marchent sur un tableau +% à coefficients formels !! +% +% ********************************************************* +% +% +% res=polyval_sym(P,x) : renvoie la valeur de P en x (en formel) : +% si n est la longueur de P : +% res=P(1)*x^(n-1)*+P(2)*x^(n-2)+...+P(n-1)*x+P(n). +% +% variables d'entrées : +% * P : polynoômes à coefficients formels +% * x : réel ou symbolique. +% variables de sortie : +% * res : valeur de P en x +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end + + +% Corps d'algorithme. +n=length(P); +if (isempty(P)) + res=0; +else + if x==0 + res=P(n); + else + res=sum(P.*x.^(n-1:-1:0)); + end +end \ No newline at end of file diff --git a/tel/TPmatlab/integration/TP3h/points_poids_gauss_diago.m b/tel/TPmatlab/integration/TP3h/points_poids_gauss_diago.m new file mode 100644 index 0000000..2b445c9 --- /dev/null +++ b/tel/TPmatlab/integration/TP3h/points_poids_gauss_diago.m @@ -0,0 +1,126 @@ +function [points,poids]=points_poids_gauss_diago(c,n) + +% points_poids_gauss_diago : calcul des poids et des points de Gauss par diagonalisation +% +% ********************************************************* +% +% [points,poids]=points_poids_gauss_diago(c,n) renvoie les points D_i et les points x_i +% des formules d'intégration à l+1 points pour 0<=l<=n, pour +% l'une des quatre méthodes : Legendre, Tchebytchev, Hermite ou Laguerre, +% calculés par diagonalisation d'une matrice tridiagonale. +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : méthode de Gauss-Legendre, +% c=2 : méthode de Gauss-Tchebychev, +% c=3 : méthode de Gauss-Hermite, +% c=4 : méthode de Gauss-Laguerre; +% +% variables de sortie +% * points : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, points(l+1,1:l+1) contient +% les l+1 points x_i de la formule (à l+1 points) +% * poids : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient +% les l+1 poids D_i de la formule (à l+1 points) +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end + +% Corps d'algorithme. +points=zeros(n+1,n+1); +poids=zeros(n+1,n+1); +switch c + case 1 + points(1,1)=0; + poids(1,1)=2; + case 2 + points(1,1)=0; + poids(1,1)=pi; + case 3 + points(1,1)=0; + poids(1,1)=sqrt(pi); + case 4 + points(1,1)=1; + poids(1,1)=1; +end +if (n>=1) + for k=1:n + if (c<=3) + B=zeros(1,k+1); + else + B=2*(0:k)+1; + end + switch c + case 1 + gamma=1./(2*sqrt(1-1./(4*(1:k).^2))); + case 2 + gamma=[1/sqrt(2) ones(1,k-1)/2]; + case 3 + gamma=sqrt((1:k)/2); + case 4 + gamma=(1:k); + end + auxi=gallery('tridiag',gamma,B,gamma); + if (k<=4) + [Q,R]=eigs(auxi,k+1); + else + [Q,R]=eig(full(auxi)); + end + R=diag(R); + U=(sqrt(sum(Q.^2))); + switch c + case 1 + auxik=2; + case 2 + auxik=pi; + case 3 + auxik=sqrt(pi); + case 4 + auxik=1; + end + Q=auxik*(Q(1,:)./U).^2; + [auxipoints,ind]=sort((flipud(R))'); + auxi=fliplr(Q); + auxipoids=auxi(ind); + if (c<=3) + p=fix(k/2); + if ~(rem(k,2)) + auxi=(auxipoints(p+2:k+1)-auxipoints(p:-1:1))/2; + auxipoints(p+2:k+1)=auxi; + auxipoints(1:p)=-fliplr(auxi); + auxipoints(p+1)=0; + auxi=(auxipoids(p+2:k+1)+auxipoids(p:-1:1))/2; + auxipoids(p+2:k+1)=auxi; + auxipoids(1:p)=fliplr(auxi); + else + auxi=(auxipoints(p+2:k+1)-auxipoints(p+1:-1:1))/2; + auxipoints(p+2:k+1)=auxi; + auxipoints(1:p+1)=-fliplr(auxi); + auxi=(auxipoids(p+2:k+1)+auxipoids(p+1:-1:1))/2; + auxipoids(p+2:k+1)=auxi; + auxipoids(1:p+1)=fliplr(auxi); + end + end + points(k+1,1:k+1)=auxipoints; + poids(k+1,1:k+1)=auxipoids; + end +end + diff --git a/tel/TPmatlab/integration/TP3h/pointspoids_diago.mat b/tel/TPmatlab/integration/TP3h/pointspoids_diago.mat new file mode 100644 index 0000000..1a5ad09 Binary files /dev/null and b/tel/TPmatlab/integration/TP3h/pointspoids_diago.mat differ diff --git a/tel/TPmatlab/integration/TP3i/demo_stocke_points_poids.m b/tel/TPmatlab/integration/TP3i/demo_stocke_points_poids.m new file mode 100644 index 0000000..4e836b3 --- /dev/null +++ b/tel/TPmatlab/integration/TP3i/demo_stocke_points_poids.m @@ -0,0 +1,84 @@ +% le script demo_stocke_points_poids permet de stocker sur le fichier pointspoids_comp +% les différents points et de poids de gauss (pour les quatre méthodes) qui sont : +% - calculés manuellement (voir TP3b) +% - calculés de façon symbolique (voir TP3g) +% - calculés par diagonalisation (voir TP3h). +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% points_poids_gauss, points_poids_gauss_symb, points_poids_gauss_diago +% +% ********************************************************* + + +clear nAmax nBmax nCmax ... + pointsleg poidsleg ... + pointslegsymb poidslegsymb ... + pointslegdiag poidslegdiag ... + pointstch poidstch ... + pointstchsymb poidstchsymb ... + pointstchdiag poidstchdiag ... + pointsher poidsher ... + pointshersymb poidshersymb ... + pointsherdiag poidsherdiag ... + pointslag poidslag ... + pointslagsymb poidslagsymb ... + pointslagdiag poidslagdiag ... + +disp('méthode directe : '); +nAmax=input('entrez l''entier maximal : ') ; +disp('méthode symbolique : '); +nBmax=input('entrez l''entier maximal : ') ; +disp('méthode par diagonalisation : '); +nCmax=input('entrez l''entier maximal : ') ; + +% calcul des poids et points de Gauss-Legendre. +% calcul manuel : +[pointsleg,poidsleg]=points_poids_gauss(1,nAmax); +% calcul symbolique : +[pointslegsymb,poidslegsymb]=points_poids_gauss_symb(1,nBmax); +% calcul par diagonalisation : +[pointslegdiag,poidslegdiag]=points_poids_gauss_diago(1,nCmax); + +% calcul des poids et points de Gauss-Tchebychev. +% calcul manuel : +[pointstch,poidstch]=points_poids_gauss(2,nAmax); +% calcul symbolique : +[pointstchsymb,poidstchsymb]=points_poids_gauss_symb(2,nBmax); +% calcul par diagonalisation : +[pointstchdiag,poidstchdiag]=points_poids_gauss_diago(2,nCmax); + +% calcul des poids et points de Gauss-Hermite. +% calcul manuel : +[pointsher,poidsher]=points_poids_gauss(3,nAmax); +% calcul symbolique : +[pointshersymb,poidshersymb]=points_poids_gauss_symb(3,nBmax); +% calcul par diagonalisation : +[pointsherdiag,poidsherdiag]=points_poids_gauss_diago(3,nCmax); + +% calcul des poids et points de Gauss-Laguerre. +% calcul manuel : +[pointslag,poidslag]=points_poids_gauss(4,nAmax); +% calcul symbolique : +[pointslagsymb,poidslagsymb]=points_poids_gauss_symb(4,nBmax); +% calcul par diagonalisation : +[pointslagdiag,poidslagdiag]=points_poids_gauss_diago(4,nCmax); + +%Stockage sur fichier +save pointspoids_comp ... + nAmax nBmax nCmax ... + pointsleg poidsleg ... + pointslegsymb poidslegsymb ... + pointslegdiag poidslegdiag ... + pointstch poidstch ... + pointstchsymb poidstchsymb ... + pointstchdiag poidstchdiag ... + pointsher poidsher ... + pointshersymb poidshersymb ... + pointsherdiag poidsherdiag ... + pointslag poidslag ... + pointslagsymb poidslagsymb ... + pointslagdiag poidslagdiag ... + + diff --git a/tel/TPmatlab/integration/TP3i/demo_trace_compar_gauss.m b/tel/TPmatlab/integration/TP3i/demo_trace_compar_gauss.m new file mode 100644 index 0000000..2c7ebce --- /dev/null +++ b/tel/TPmatlab/integration/TP3i/demo_trace_compar_gauss.m @@ -0,0 +1,134 @@ +% le script trace_compar_gauss permet de comparer les différents +% calculs de points et de poids de gauss en comparant l'intégrale exacte +% et l'intégrale approchée d'une fonction selon que les poids et points +% ont été calculés : +% * par méthode directe; +% * par méthode symbolique; +% * ou par méthode par diagonalisation, +% et ce pour chacune des quatre méthodes définies par c. +% les poids et les points seront préalablements calculés selon chacune des trois méthodes +% dans les tableaux : +% pointsleg poidsleg ... +% pointslegsymb poidslegsymb ... +% pointslegdiag poidslegdiag ... +% pointstch poidstch ... +% pointstchsymb poidstchsymb ... +% pointstchdiag poidstchdiag ... +% pointsher poidsher ... +% pointshersymb poidssherymb ... +% pointsherdiag poidsherdiag ... +% pointslag poidslag ... +% pointslagsymb poidslagsymb ... +% pointslagdiag poidslagdiag ... +% voir script demo_stocke_points_poids +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% saisiefonction, int_gauss +% +% ********************************************************* + + + + +% logarithme de l'erreur nulle prise égale à log10zero +log10zero=-35; + + +clear xA yA xB yB xC yC Iexa erA erB erC; +disp('choix de la méthode : ') +c=input('entrez 1, 2, 3 ou 4 : '); +disp(' '); +disp('Entrée de la fonction r telle que r=f/w sur (a,b) où w est le poids ') ; +disp(' '); +fcn=saisiefonction; +Iexa=input('entrez la valeur exacte de l''intégrale de r=f/w sur (a,b) : ') ; +erA=zeros(1,nAmax+1); +yA=zeros(1,nAmax+1); +erB=zeros(1,nBmax+1); +yB=zeros(1,nBmax+1); +erC=zeros(1,nCmax+1); +yC=zeros(1,nCmax+1); +switch c + case 1 + pointsA=pointsleg; + poidsA=poidsleg; + pointsB=pointslegsymb; + poidsB=poidslegsymb; + pointsC=pointslegdiag; + poidsC=poidslegdiag; + case 2 + pointsA=pointstch; + poidsA=poidstch; + pointsB=pointstchsymb; + poidsB=poidstchsymb; + pointsC=pointstchdiag; + poidsC=poidstchdiag; + case 3 + pointsA=pointsher; + poidsA=poidsher; + pointsB=pointshersymb; + poidsB=poidshersymb; + pointsC=pointsherdiag; + poidsC=poidsherdiag; + case 4 + pointsA=pointslag; + poidsA=poidslag; + pointsB=pointslagsymb; + poidsB=poidslagsymb; + pointsC=pointslagdiag; + poidsC=poidslagdiag; +end + + +for k=0:nAmax + auxi=int_gauss(k,pointsA,poidsA,fcn); + erA(k+1)=abs(auxi-Iexa); +end +for k=0:nBmax + auxi=int_gauss(k,pointsB,poidsB,fcn); + auxi=vpa(auxi)-Iexa; + if ~(isnumeric(auxi)) + auxi=double(auxi); + end + erB(k+1)=abs(auxi); +end +for k=0:nCmax + auxi=int_gauss(k,pointsC,poidsC,fcn); + erC(k+1)=abs(auxi-Iexa); +end +indnul=find(erA==0); +indnonnul=find(erA~=0); +yA(indnul)=log10zero; +yA(indnonnul)=log10(erA(indnonnul)); +indnul=find(erB==0); +indnonnul=find(erB~=0); +yB(indnul)=log10zero; +yB(indnonnul)=log10(erB(indnonnul)); +indnul=find(erC==0); +indnonnul=find(erC~=0); +yC(indnul)=log10zero; +yC(indnonnul)=log10(erC(indnonnul)); +xA=log10(1:nAmax+1); +xB=log10(1:nBmax+1); +xC=log10(1:nCmax+1); +clf; +hold on; +plot(xA,yA,'b',xB,yB,'r',xC,yC,'g'); +xlabel('logarithme 10 du nombre de points d''intégration'); +ylabel('logarithme 10 de l''erreur d''intégration'); +legend('directe','symbolique','diagonalisation',0); +hold off; + + + + + + + + + + + + diff --git a/tel/TPmatlab/integration/TP3i/verifie_nullite_gauss.m b/tel/TPmatlab/integration/TP3i/verifie_nullite_gauss.m new file mode 100644 index 0000000..8288744 --- /dev/null +++ b/tel/TPmatlab/integration/TP3i/verifie_nullite_gauss.m @@ -0,0 +1,103 @@ +function [eta1,eta2]=verifie_nullite_gauss(c,n) + + +% verifie_nullite_gauss : calcul de l'orthogonalité et de la normalisation des polynômes calculés numériquement +% +% ********************************************************* +% +% [eta1,eta2]=verifie_nullite_gauss(c,n) : +% vérifie numériquement que les polynômes de Gauss P_0, ..., P_n +% sont orthogonaux avec la condition de normalité : +% P_r(1)=1 pour Legendre et Tchebytchev +% P_r(0)=1 pour Laguerre +% le coefficient domimant de P_r est 2^{r} pour Hermite +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : Gauss-Legendre, +% c=2 : Gauss-Tchebytchev, +% c=3 : Gauss-Hermite, +% c=4 : Gauss-Laguerre. +% variables de sortie +% * eta1= max_{0<=r<=n} |P_r(1)-1| pour Legendre et Tchebychev +% eta1= max_{0<=r<=n} |P_r(0)-1| pour Laguerre +% eta1= max_{0<=r<=n} |2^r-coefficient dominant(P_r)|=0 pour Hermite +% et +% eta2=max_{0<=r,s<=n, r~=s}|| +% (calculé par intégration en symbolique) +% +% eta1 et eta2 doivent être nuls. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% calcul_poly_gauss +% +% ********************************************************* + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end +if ~((c==1) | (c==2) | (c==3) | (c==4)) + error('c doit être un entier égal à 1, 2, 3 ou 4'); +end + +% Corps d'algorithme. +P=calcul_poly_gauss(c,n); +ermax=0; +if (c==1) | (c==2) + for r=0:n + er=abs(polyval(P(r+1,1:r+1),1)-1); + ermax=max(er,ermax); + end +end +if (c==3) + for r=0:n + er=abs(P(r+1,1)-2^r); + ermax=max(er,ermax); + end +end +if (c==4) + for r=0:n + er=abs(polyval(P(r+1,1:r+1),0)-1); + ermax=max(er,ermax); + end +end +eta1=ermax; +syms x a b; +switch c + case 1 + a=sym(-1); + b=sym(1); + w=sym(1); + case 2 + a=sym(-1); + b=sym(1); + w=1/sqrt(1-x^2); + case 3 + a=-sym(inf); + b=sym(inf); + w=exp(-x^2); + case 4 + a=sym(0); + b=sym(+inf); + w=exp(-x); +end +ermax=0; +for r=0:n + for s=0:r-1 + R=conv(P(r+1,1:r+1),P(s+1,1:s+1)); + f=poly2sym(R); + er=double(abs(simplify(int(f*w,a,b)))); + ermax=max(er,ermax); + end +end +eta2=ermax; diff --git a/tel/TPmatlab/integration/TP3i/verifie_nullite_gauss_symb.m b/tel/TPmatlab/integration/TP3i/verifie_nullite_gauss_symb.m new file mode 100644 index 0000000..f11e97f --- /dev/null +++ b/tel/TPmatlab/integration/TP3i/verifie_nullite_gauss_symb.m @@ -0,0 +1,120 @@ +function [eta1,eta2]=verifie_nullite_gauss_symb(c,n) + + +% verifie_nullite_gauss_symb : calcul de l'orthogonalité et de la normalisation des polynômes en symbolique +% +% ********************************************************* +% +% [eta1,eta2]=verifie_nullite_gauss_symb(c,n) : +% vérifie en symbolique que les polynômes de Gauss P_0, ..., P_n +% sont orthogonaux avec la condition de normalité : +% P_r(1)=1 pour Legendre et Tchebytchev +% P_r(0)=1 pour Laguerre +% le coefficient domimant de P_r est 2^{r} pour Hermite +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : Gauss-Legendre, +% c=2 : Gauss-Tchebytchev, +% c=3 : Gauss-Hermite, +% c=4 : Gauss-Laguerre. +% variables de sortie eta1 et eta2 où en symbolique, +% eta1=0 si max_{0<=r<=n} |P_r(1)-1|=0 pour Legendre et Tchebychev +% si max_{0<=r<=n} |P_r(0)-1|=0 pour Laguerre +% si max_{0<=r<=n} |2^r-coefficient dominant(P_r)|=0 pour Hermite +% et +% eta2=0 si max_{0<=r,s<=n, r~=s}||=0 +% (en symbolique) +% +% eta1 et eta2 doivent être nuls. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% calcul_poly_gauss_symb, polyval_symb +% +% ********************************************************* + + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if (fix(n)~=n) | (n<0) + error('l''entier n doit être un entier naturel'); +end +if ~((c==1) | (c==2) | (c==3) | (c==4)) + error('c doit être un entier égal à 1, 2, 3 ou 4'); +end + +% Corps d'algorithme. +P=calcul_poly_gauss_symb(c,n); +ind=0; +if (c==1) | (c==2) + for r=0:n + auxi=polyval_symb(P(r+1,1:r+1),1)-sym(1); + indloc=~(auxi==0); + ind=max(ind,indloc); + end +end +if (c==3) + for r=0:n + auxi=P(r+1,1); + auxi=auxi-2^(sym(r)); + indloc=~(auxi==0); + ind=max(ind,indloc); + end +end +if (c==4) + for r=0:n + auxi=polyval_symb(P(r+1,1:r+1),0)-sym(1); + indloc=~(auxi==0); + ind=max(ind,indloc); + end +end +eta1=ind; +disp('vérification initiale'); +if ~(ind) + disp('ok !!'); +else + disp('erreur !!'); +end +syms x a b; +switch c + case 1 + a=sym(-1); + b=sym(1); + w=sym(1); + case 2 + a=sym(-1); + b=sym(1); + w=1/sqrt(1-x^2); + case 3 + a=-sym(inf); + b=sym(inf); + w=exp(-x^2); + case 4 + a=sym(0); + b=sym(+inf); + w=exp(-x); + end +ind=0; +for r=0:n + for s=0:r-1 + f=polyval_symb(P(r+1,1:r+1),x); + g=polyval_symb(P(s+1,1:s+1),x); + auxi=simplify(int(f*g*w,a,b)); + indloc=~(auxi==0); + ind=max(indloc,ind); + end + disp(['vérification pour r=',int2str(r)]); + if ~(ind) + disp('ok !!'); + else + disp('erreur !!'); + end +end +eta2=ind; diff --git a/tel/TPmatlab/integration/TP3i/verifie_ordre_gauss.m b/tel/TPmatlab/integration/TP3i/verifie_ordre_gauss.m new file mode 100644 index 0000000..d3a5ecd --- /dev/null +++ b/tel/TPmatlab/integration/TP3i/verifie_ordre_gauss.m @@ -0,0 +1,63 @@ +function eta=verifie_ordre_gauss(c,n,points,poids) + + +% verifie_ordre_gauss : vérification de l'ordre des formules Gaussiennes +% +% ********************************************************* +% +% eta=verifie_ordre_gauss(c,n,points,poids) calcule numériquement le max +% de |I_{n,k}-int_a ^b x^k w(x)dx| pour 0<=k<=2n+1, où I_{n,k} +% est l'intégrale approchée de x^kw(x) sur l'intervalle conventionnel (a,b) +% +% variables d'entrées : +% * n : entier naturel. +% * c définit la méthode : +% c=1 : Gauss-Legendre, +% c=2 : Gauss-Tchebytchev, +% c=3 : Gauss-Hermite, +% c=4 : Gauss-Laguerre. +% * points : tableau de taille (q+1)*(q+1) (avec n<=q) tel que +% pour tout l dans {0,...,q}, points(l+1,1:l+1) contient +% les l+1 points x_i de la formule (à l+1 points) +% * poids : tableau de taille (n+1)*(n+1) tel que +% pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient +% les l+1 poids D_i de la formule (à l+1 points) +% variable de sortie +% * eta= max_{0<=k<=2n+1} |I_{n,k}-int_a ^b x^k w(x)dx| +% +% eta doit être nul. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% expression_part_symb +% +% ********************************************************* + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=4 + error('nombre d''arguments de la fonction incorrect'); +end +% autres tests éventuels +if ~((c==1) | (c==2) | (c==3) | (c==4)) + error('c doit être un entier égal à 1, 2, 3 ou 4'); +end +p=length(poids); +if (2*n+1>=p) + error('n est trop grand'); +end + + +% Corps d'algorithme. +er=0; +for k=0:2*n+1 + I=sum(poids(k+1,1:k+1).*((points(k+1,1:k+1)).^k)); + erloc=expression_part_symb(c,k)-I; + if ~(isnumeric(erloc)) + erloc=double(vpa(erloc)); + end + erloc=abs(erloc); + er=max(er,erloc); +end +eta=er; \ No newline at end of file diff --git a/tel/TPmatlab/interpolation/TP2a/diff_div_dist.m b/tel/TPmatlab/interpolation/TP2a/diff_div_dist.m new file mode 100644 index 0000000..b3b76b5 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2a/diff_div_dist.m @@ -0,0 +1,60 @@ +function D=diff_div_dist(X,Y) +% +% +% diff_div_dist : calcul des différences divisées pour des points de supports distincts +% +% ********************************************************* +% +% D=diff_div_dist(X,Y) +% les différences divisées sont calculées par l'algorithme +% pyramidal donné en cours. +% +% variables d'entrées : +% * X : contient les valeurs x_i, pour 0 <=i<=n (deux à deux distinctes) +% * Y : contient les valeurs f(x_i), pour 0 <=i<=n +% +% variables de sortie +% * D : contient les différences divisées +% f[x_0], f[x_0,x_1], ...., f[x_0,x_1,...,x_n] +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% taille des deux variables +n=length(X)-1; +np=length(Y)-1; +if n~=np + error('les deux tableaux n''ont pas la même taille'); +end +% vérification des éléments du support deux à deux disjoints +% (si X non formel) +if isnumeric(X) + Z=sort(X); + if min(diff(Z))==0 + error('deux points du support sont égaux'); + end +end + + +% Corps d'algorithme +D=Y; +for i=1:n; + for j=n:-1:i; + D(j+1)=(D(j+1)-D(j))/(X(j+1)-X(j-i+1)); + end +end + + + diff --git a/tel/TPmatlab/interpolation/TP2a/eval_horner.m b/tel/TPmatlab/interpolation/TP2a/eval_horner.m new file mode 100644 index 0000000..7a363e7 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2a/eval_horner.m @@ -0,0 +1,49 @@ +function Z=eval_horner(T,C,D) +% +% +% eval_horner : évaluation vectorielle d'un polynôme à partir de sa forme de Newton. +% +% ********************************************************* +% +% Z=eval_horner(T,C,D) +% L'évaluation est faite par l'algorithme de Hörner donné en cours, +% Elle permet de calculer : +% p_n(T)=d_0+d_1(T-c_1)+d_2(T-c_1)(T-c_2)+...+d_n(T-c_1)...(T-c_n) +% +% variables d'entrées : +% * T : contient le tableau rectangulaire des réels où est évalué p_n +% * C : contient les centres c_i, pour 1 <=i<=n +% * D : contient les coefficients d_i, pour 0<=i<=n +% +% variables de sortie +% * Z : contient le tableau rectangulaire des images de T par p_n +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=3 + error('nombre d''arguments de la fonction incorrect'); +end +% taille des variables C et D +n=length(C); +np=length(D)-1; +if n~=np + error('les deux tableaux n''ont pas la taille adéquate'); +end + + +% Corps d'algorithme +O=ones(size(T)); +Auxi=D(n+1).*O; +for i=n-1:-1:0 + Auxi=(T-C(i+1)).*Auxi+D(i+1); +end +Z=Auxi; diff --git a/tel/TPmatlab/interpolation/TP2a/eval_horner_sca.m b/tel/TPmatlab/interpolation/TP2a/eval_horner_sca.m new file mode 100644 index 0000000..6f0dd30 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2a/eval_horner_sca.m @@ -0,0 +1,51 @@ +function z=eval_horner_sca(t,C,D) +% +% +% eval_horner_sca : évaluation scalaire d'un polynôme à partir de sa forme de Newton. +% +% ********************************************************* +% +% z=eval_horner_sca(t,C,D) +% L'évaluation est faite par l'algorithme de Hörner donné en cours, +% Elle permet de calculer : +% p_n(t)=d_0+d_1(t-c_1)+d_2(t-c_1)(t-c_2)+...+d_n(t-c_1)...(t-c_n) +% +% variables d'entrées : +% * t : contient le réel où est évalué p_n +% * C : contient les centres c_i, pour 1 <=i<=n +% * D : contient les coefficients d_i, pour 0<=i<=n +% +% variables de sortie +% * z : contient p_n(t) +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=3 + error('nombre d''arguments de la fonction incorrect'); +end +% taille des variables C et D +n=length(C); +np=length(D)-1; +if n~=np + error('les deux tableaux n''ont pas la taille adéquate'); +end + + +% Corps d'algorithme +auxi=D(n+1); +for i=n-1:-1:0 + auxi=(t-C(i+1))*auxi+D(i+1); +end +z=auxi; + + + diff --git a/tel/TPmatlab/interpolation/TP2a/test_diff_div_dist.m b/tel/TPmatlab/interpolation/TP2a/test_diff_div_dist.m new file mode 100644 index 0000000..984f048 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2a/test_diff_div_dist.m @@ -0,0 +1,32 @@ +function R=test_diff_div_dist(T,X,Y) + +% test_diff_div_dist : fonction test pour le calcul du polynôme interpolateur (vectoriel) +% +% +% ********************************************************* +% +% R=test_diff_div_dist(T,X,Y) +% calcul de la valeur du polynôme d'interpolation défini par le nuage de point X,Y +% sur le tableau T (utilisation du calcul des différences divisées et +% de l'algorithme d'évaluation de Horner). +% +% variables d'entrées : +% * T : contient le tableau rectangulaire des réels où est évalué p_n +% * X : contient les valeurs x_i, pour 0 <=i<=n (deux à deux distinctes) +% * Y : contient les valeurs f(x_i), pour 0 <=i<=n +% +% variables de sortie +% * R : contient le tableau rectangulaire des images de T par p_n, +% polynôme d'interpolation de f sur le support x_i, pour 0 <=i<=n +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div_dist, eval_horner +% +% ********************************************************* + + +D=diff_div_dist(X,Y); +R=eval_horner(T,X(1:end-1),D); + diff --git a/tel/TPmatlab/interpolation/TP2a/test_diff_div_dist_sca.m b/tel/TPmatlab/interpolation/TP2a/test_diff_div_dist_sca.m new file mode 100644 index 0000000..f612c4a --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2a/test_diff_div_dist_sca.m @@ -0,0 +1,32 @@ +function r=test_diff_div_dist_sca(t,X,Y) + +% test_diff_div_dist_sca : fonction test pour le calcul du polynôme interpolateur (scalaire) +% +% +% ********************************************************* +% +% R=test_diff_div_dist_sca(t,X,Y) +% calcul de la valeur du polynôme d'interpolation défini par le nuage de point X,Y +% en t (utilisation du calcul des différences divisées et +% de l'algorithme d'évaluation de Hörner). +% +% variables d'entrées : +% * t : contient le réel où est évalué p_n +% * X : contient les valeurs x_i, pour 0 <=i<=n (deux à deux distinctes) +% * Y : contient les valeurs f(x_i), pour 0 <=i<=n +% +% variables de sortie +% * r : contient p_n(t)où p_n est le +% polynôme d'interpolation de f sur le support x_i, pour 0 <=i<=n +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div_dist, eval_horner_sca +% +% ********************************************************* + + +D=diff_div_dist(X,Y); +r=eval_horner_sca(t,X(1:end-1),D); + diff --git a/tel/TPmatlab/interpolation/TP2b/lect_courb.m b/tel/TPmatlab/interpolation/TP2b/lect_courb.m new file mode 100644 index 0000000..b991130 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2b/lect_courb.m @@ -0,0 +1,59 @@ +function lect_courb(nb_chif) +% lit les coordonnées des points d'une courbe +% +% variables d'entrée et fonctionnement +% Tant que l'utilisateur clique gauche, +% le système affiche les coordonnées du point +% choisi à la souris. +% nb_chif désigne le nombre de décimales demandées, +% afin de ne pas encombrer l'écran en cas de lectures +% multiples. On impose 0<=nb_chif<=3, aisément modifiable. +% On sort en cliquant droit. +% +% variables de sortie +% affichage des coordonnées du point cliqué. +% +% +% ************ Fonctions connexes utilisées *************** +% +% ********************************************************* + +% Contrôles d'entrée +% nombre d'arguments +if nargin~=1 + error('Nombre d''argument incorrect. Revoir le source.'); +end + + +% autres tests éventuels +if (nb_chif<0)|(nb_chif>3) + error('Champ incorrect.Modifier le source si besoin'); +end + + +% Corps d'algorithme + +% traitement proprement dit +bout=1; +while bout==1 + [x,y,bout]=ginput(1); + % réglage des sorties + entx=length(num2str(floor(x))); + enty=length(num2str(floor(y))); + fx=['%' num2str(entx) '.' num2str(nb_chif) 'f']; + fy=['%' num2str(enty) '.' num2str(nb_chif) 'f']; + text(x,y,['x=' num2str(x,fx) ' y=' num2str(y,fy)]); +end + + +% variante pour la sortie des résultats +% L'utilisateur pourra afficher seulement le numéro des points +% saisis sur la courbe; il remplira en meme temps une table +% des coordonnées qu'il peut ajouter comme paramètre de sortie +% de cette nouvelle fonction. + + +% fin de fonction() + + + diff --git a/tel/TPmatlab/interpolation/TP2c/etude_support.m b/tel/TPmatlab/interpolation/TP2c/etude_support.m new file mode 100644 index 0000000..560dcc3 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2c/etude_support.m @@ -0,0 +1,110 @@ +function etude_support(a,b,n) +% compare l'effet sur la qualité de l'interpolation +% des points de Tchebychev et +% de points régulièrement répartis dans [a,b]. +% +% variables d'entrée +% n est un entier naturel +% tel que x0,...,xn sont les points d'interpolation +% choisis dans [a,b] selon les deux variantes à comparer. +% +% variables de sortie +% On visualise; on pourra capter l'image produite +% et éventuellement renvoyer son handle. +% +% +% ************ Fonctions auxiliaires utilisées *************** +% +% diff_div_dist, eval_horner +% +% ********************************************************* + +% Paramètres modifiables +format long e; +nmax=19; % nb de points d'interpolation +N=100; % détermine la taille des vecteurs à représenter +pas=(b-a)/N; +% Contrôles d'entrée + +% nombre d'arguments +if nargin~=3 + error('Nombre des arguments incorrects'); +end + +% autres tests éventuels +if n>nmax + error('Modifiez le paramètre nmax en tete de fichier') +end +if b<=a + error('Choisissez un intervalle conventionnel') +end + +% Corps d'algorithme + +% préparation des deux familles de points de support +i=0:n; +% points régulièrement répartis +x1=a+(b-a)/n.*i; +% points de Tchebychev +x2=(b-a)/2*cos(pi/2*(2.*i+1)/(n+1))+(b+a)/2; + +% détermination des images d'un vecteur x +% par la fonction proposée; +% on pourrait rajouter un champ de saisie +% de la fonction à interpoler. +X=a:pas:b; +Y=exp(X); + +% détermination des polynomes d'interpolation +% et des images du vecteur x obtenues par eux. +y1=exp(x1);y2=exp(x2); +d1=diff_div_dist(x1,y1);d2=diff_div_dist(x2,y2); +Y1=eval_horner(X,x1(1:n),d1); +Y2=eval_horner(X,x2(1:n),d2); + +% calcul des écarts entre fonction exacte +% et polynomes d'interpolation +E1=abs(Y-Y1);E2=abs(Y-Y2); + +% calibrage des images +maxy=max(max(E1),max(E2)); +bordx=2*(b-a)/10; +lim=[a-bordx b+bordx 0 1.1*maxy]; + +clf; +% représentation de E1 et E2 +subplot(3,1,1) +plot(X,E1,'b:'); +axis(lim); +title(['Erreurs d''interpolation à ' num2str(n+1) ' points pour exponentielle']); +legend('Cas de points régulièrement répartis',0); + +subplot(3,1,2) +plot(X,E2,'r-'); +axis(lim); +legend('Cas de points de Tchebychev',0); + +subplot(3,1,3) +plot(X,E1,'b:'); +axis(lim); +legend('Erreurs comparées',0); +hold on; +plot(X,E2,'r-'); +axis(lim); + + + + + + + + + + + +% sortie des résultats + + + +% fin de fonction() + diff --git a/tel/TPmatlab/interpolation/TP2d/demo_test_compar_methode.m b/tel/TPmatlab/interpolation/TP2d/demo_test_compar_methode.m new file mode 100644 index 0000000..37ce123 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2d/demo_test_compar_methode.m @@ -0,0 +1,90 @@ +% script demo_test_compar_methode : comparaison des différentes méthodes pour calculer un polynôme d'interpolation +% +% +% Ce script permet de comparer, en terme d'évaluation numérique et graphique +% le polynôme interpolateur p_n d'une fonction f sur le support (x_i) où, pour 0<=i<=n, +% on a x_i=a+i*(b-a)/n, avec a, b, n et f donnée, en fonction des trois méthodes : +% a) calcul en utilisant la forme de Newton et les différences divisées, +% b) calcul à partir des polynômes de Lagranges l_i. +% c) calcul à partir de la forme canonique de p_n +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% saisiefonction, diff_div_dist, eval_horner, +% newton_cano, poly_eval_lagr +% +% ********************************************************* + +clear all; +a=input('Entrez l''abscisse minimale a : '); +b=input('Entrez l''abscisse maximale b : '); +pgra=input('Entrez le nombre de points pour le graphique : '); +f=saisiefonction; +n=input('Entrez l''entier n positif (égal au nombre de points oté de un) : '); +beton=0; +choixpoint='o'; +xgra=a:(b-a)/pgra:b; +ygra=feval(f,xgra); +while(n>=0) + if strcmp(choixpoint,'o') + choixpoint=input('Afficher le support ? entrez o ou n : ','s'); + end + switch beton + case 0 + choix=input('Conserver le graphique par calcul à partir de la forme canonique ? entrez o ou n : ','s'); + if strcmp(choix,'n') + beton=1; + end + case 1 + choix=input('Conserver le graphique par calcul à partir des polynômes de Lagrange ? entrez o ou n : ','s'); + if strcmp(choix,'n') + beton=2; + end + end + if (n==0) + X=(a+b)/2; + else + X=a:(b-a)/n:b; + end + Y=feval(f,X); + D=diff_div_dist(X,Y); + ygranewton=eval_horner(xgra,X(1:end-1),D); + ygralagrange=poly_eval_lagr(xgra,X,Y); + P=newton_cano(D,X(1:end-1)); + ygracanonique=polyval(P,xgra); + clf; + hold on; + if (strcmp(choixpoint,'o')) + switch beton + case 0 + plot(X,Y,'go',xgra,ygra,xgra,ygranewton,xgra,ygralagrange,xgra,ygracanonique); + legend('Points du support','Fonction exacte','Newton','Lagrange','Canonique',0); + case 1 + plot(X,Y,'go',xgra,ygra,xgra,ygranewton,xgra,ygralagrange); + legend('Points du support','Fonction exacte','Newton','Lagrange',0); + case 2 + plot(X,Y,'go',xgra,ygra,xgra,ygranewton); + legend('Points du support','Fonction exacte','Newton',0); + end + else + switch beton + case 0 + plot(xgra,ygra,xgra,ygranewton,xgra,ygralagrange,xgra,ygracanonique); + legend('Fonction exacte','Newton','Lagrange','Canonique',0); + case 1 + plot(xgra,ygra,xgra,ygranewton,xgra,ygralagrange); + legend('Fonction exacte','Newton','Lagrange',0); + case 2 + plot(xgra,ygra,xgra,ygranewton); + legend('Fonction exacte','Newton',0); + end + end + if exist(f)==1 + title(['interpolation de la fonction ',formula(f),' avec ',int2str(n+1),' points']); + else + title(['interpolation de la fonction ',f,' avec ',int2str(n+1),' points']); + end + hold off; + n=input('pour continuer, entrez un nombre n positif (pour arrêter, le choisir négatif) : '); +end \ No newline at end of file diff --git a/tel/TPmatlab/interpolation/TP2d/newton_cano.m b/tel/TPmatlab/interpolation/TP2d/newton_cano.m new file mode 100644 index 0000000..81055ee --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2d/newton_cano.m @@ -0,0 +1,54 @@ +function B=newton_cano(A,C) +% +% +% newton_cano : passage de la forme de newton à la forme canonique. +% +% ********************************************************* +% +% B=newton_cano(A,C) +% L'algorithme repose sur l'algorithme d'Hörner réitéré n fois, avec +% décalage successif des centres. Cette fonction permet de passer de +% la forme de Newton : +% p_n(x)=a_0+a_1(x-c_1)+a_2(x-c_1)(x-c_2)+...+a_n(x-c_1)...(x-c_n) +% à la forme canonique +% p_n(x)=b_0+b_1x+b_2x^2+...+b_nx^n +% +% variables d'entrées : +% * C : contient les centres c_i, pour 1 <=i<=n +% * A : contient les coefficients a_i, pour 0<=i<=n +% +% variables de sortie +% * B : contient les coefficients b_n, b_(n-1),...,b_0 +% (dans cet ordre là) +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% taille des variables X et D +n=length(C); +np=length(A)-1; +if n~=np + error('les deux tableaux n''ont pas la taille adéquate'); +end + + +% Corps d'algorithme +Auxi=A; +for j=1:n + for i=n-1:-1:j-1 + Auxi(i+1)=Auxi(i+1)-C(i-j+2)*Auxi(i+2); + end +end +B=fliplr(Auxi); diff --git a/tel/TPmatlab/interpolation/TP2d/poly_eval_lagr.m b/tel/TPmatlab/interpolation/TP2d/poly_eval_lagr.m new file mode 100644 index 0000000..718abf6 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2d/poly_eval_lagr.m @@ -0,0 +1,50 @@ +function Z=poly_eval_lagr(T,X,Y) + +% poly_lagrange : calcul du polynôme de Lagrange en T (vectoriel) +% +% ********************************************************* +% +% Z=poly_eval_lagr(T,X,Y) +% calcul de la valeur du polynôme p_n d'interpolation défini par le nuage de point X,Y +% en T (à partir des polynômes de Lagrange l_i). +% +% variables d'entrées : +% * T : contient le tableau ligne des réels où on évalue p_n +% * X : contient les centres x_i, pour 0<=i<=n (deux à deux distincts) +% * Y : contient les valeurs f(x_i), pour 0<=i<=n +% +% variables de sortie +% * Z : contient l'image du tableau T par p_n (tableau ligne) +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% poly_lagrange +% +% ********************************************************* +% + + + + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=3 + error('nombre d''arguments de la fonction incorrect'); +end +na=length(X); +nb=length(Y); +if (na~=nb) + error('X et Y n''ont pas la même taille'); +end +n=na-1; + +% Corps d'algorithme +Auxi=Y(1)*poly_lagrange(0,T,X); +for i=2:n+1 + Auxi=Auxi+Y(i)*poly_lagrange(i-1,T,X); +end +Z=Auxi; diff --git a/tel/TPmatlab/interpolation/TP2d/poly_lagrange.m b/tel/TPmatlab/interpolation/TP2d/poly_lagrange.m new file mode 100644 index 0000000..66f5871 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2d/poly_lagrange.m @@ -0,0 +1,53 @@ +function Z=poly_lagrange(i,T,X) + +% poly_lagrange : calcul du i-ième polynôme de Lagrange en T (vectoriel) +% +% ********************************************************* +% +% z=poly_lagrange_sca(i,T,X) +% Calcul de l_i(t) défini par +% l_i(T)=((T-x_0)/(x_i-x_0))...((T-x_(i-1))/(x_i-x_(i-1))) +% ((T-x_(i+1))/(x_i-x_(i+1)))... +% ((T-x_n)/(x_i-x_n)) +% +% variables d'entrées : +% * i : contient l'indice i, compris entre 0 et n +% * T : contient le tableau ligne des réels où on évalue l_i +% * X : contient les centres x_i, pour 0<=i<=n (deux à deux distincts) +% +% variables de sortie +% * Z : contient l'image du tableau T par l_i (tableau ligne) +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=3 + error('nombre d''arguments de la fonction incorrect'); +end +n=length(X)-1; +if (i<0) | (i>n) + error('i n''appartient pas à l''ensemble {0,...,n}') +end +% vérification des éléments du support deux à deux disjoints +Z=sort(X); +if min(diff(Z))==0 + error('deux points du support sont égaux'); +end + +% Corps d'algorithme +q=length(T); +y=X(1:i); +y(i+1:n)=X(i+2:n+1); +Auxi=X(i+1)-X(1:i); +Auxi(i+1:n)=X(i+1)-X(i+2:n+1); +B=ones(n,1)*T-y'*ones(1,q); +Z=prod(((1./Auxi)'*ones(1,q)).*B,1); diff --git a/tel/TPmatlab/interpolation/TP2d/poly_lagrange_sca.m b/tel/TPmatlab/interpolation/TP2d/poly_lagrange_sca.m new file mode 100644 index 0000000..b0fe2e2 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2d/poly_lagrange_sca.m @@ -0,0 +1,51 @@ +function z=poly_lagrange_sca(i,t,X) + +% poly_lagrange_sca : calcul du i-ième polynôme de Lagrange en t (scalaire) +% +% ********************************************************* +% +% z=poly_lagrange_sca(i,t,X) +% Calcul de l_i(t) défini par +% l_i(t)=((t-x_0)/(x_i-x_0))...((t-x_(i-1))/(x_i-x_(i-1))) +% ((x-x_(i+1))/(x_i-x_(i+1)))... +% ((x-x_n)/(x_i-x_n)) +% +% variables d'entrées : +% * i : contient l'indice i, compris entre 0 et n +% * t : contient le réel t où on évalue l_i +% * X : contient les centres x_i, pour 0<=i<=n (deux à deux distincts) +% +% variables de sortie +% * z : contient l_i(t) +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=3 + error('nombre d''arguments de la fonction incorrect'); +end +n=length(X)-1; +if (i<0) | (i>n) + error('i n''appartient pas à l''ensemble {0,...,n}') +end +% vérification des éléments du support deux à deux disjoints +Z=sort(X); +if min(diff(Z))==0 + error('deux points du support sont égaux'); +end + +% Corps d'algorithme +y=X(1:i); +y(i+1:n)=X(i+2:n+1); +Auxi=X(i+1)-X(1:i); +Auxi(i+1:n)=X(i+1)-X(i+2:n+1); +z=prod((t-y)./Auxi); diff --git a/tel/TPmatlab/interpolation/TP2e/elem_simp.m b/tel/TPmatlab/interpolation/TP2e/elem_simp.m new file mode 100644 index 0000000..9786230 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2e/elem_simp.m @@ -0,0 +1,83 @@ +function [ch_pol,ch_elem,coef]=elem_simp(x,N) +% décompose certaines fractions rationnelles en éléments simples +% +% variables d'entrée +% x est le vecteur réel des poles de la fraction rationnelle, +% dont sera issu par la fonction poly, le dénominateur D(x). +% N est le numérateur de la fraction rationnelle, +% passé comme un polynome matlab, c'est-à-dire un vecteur. +% +% variables de sortie +% ch_pol est une chaine représentant la partie polynomiale, quotient +% obtenu dans la division euclidienne de N par D(x)=(x-x1)...(x-xn). +% ch_elem est une chaine représentant les éléments simples +% au sens ordinaire en fonction de coefficients fournis dans coef. +% coef peut etre forcé à une écriture symbolique (ici rationnelle). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Fonctions connexes appelées +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Tests d'entrée généraux +if nargin~=2 + error('Nombre de chapms incorrect'); +end +if isreal(x)==0 + error('Les poles doivent etre réels pour l''instant'); +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Mise sous forme conventionnelle de la fraction +% calcul du dénominateur +D=poly(x); +[Q,R]=deconv(N,D); +% désormais on travaille avec R/D qui correspond à l'étude +% proposée au début dans le tp. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% Premier élément du résultat +ch_pol=poly2str(Q,'x'); + + +% recherche de l'existence de poles multiples dans x +% tri de x pour faciliter le repérage des redondances +x=sort(x); +% détection du nombre de poles multiples +pol_mult=sum(diff(x)==0); + +% ALGORITHME GENERAL + +switch pol_mult +case 0, + % il s'agit de l'étude demandée dans ce tp. + + der=polyder(D); + coef=polyval(R,x)./polyval(der,x); + ch_elem=[]; + for i=1:size(x,2) + switch sign(x(i)) + case 1, + ch=['/(x-' num2str(x(i)) ')']; + case -1, + ch=['/(x+' num2str(-x(i)) ')']; + case 0, + ch='/x'; + end + ch_elem=[ch_elem ' + a' num2str(i) ch]; + end + % à supprimer éventuellement par transformation en commentaire! + coef=sym(coef,'r'); + +otherwise + % étude à compléter par l'utilisateur... + % Voir indications dans le fichier de documentation. + disp('Il existe au moins un pole multiple. Développez l''étude générale.') +end + + + + + + diff --git a/tel/TPmatlab/interpolation/TP2f/demo_limite_diff_div.m b/tel/TPmatlab/interpolation/TP2f/demo_limite_diff_div.m new file mode 100644 index 0000000..a106c43 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/demo_limite_diff_div.m @@ -0,0 +1,121 @@ +% script demo_limite_diff_div: calcul et comparaisons de différences divisées par limite formelle et dérivation +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div +% +% ********************************************************* +% +% + +% corps d'algorithme +clear all; +syms x0 x1 x2; +syms y0 y1 y2; +syms x; +f=sym('f(x)'); +fp=diff(f); +fpp=diff(f,2); +y0=subs(f,x,x0); +y1=subs(f,x,x1); +y2=subs(f,x,x2); + +% calcul avec les trois valeurs distinctes +disp(' '); +disp(' '); +disp('calcul avec les trois valeurs distinctes'); +disp('valeurs des trois différences divisées f[x0],f[x0,x1],f[x0,x1,x2]'); +disp(' '); +D=diff_div([x0 x1 x2],[y0 y1 y2]); +pretty(D); +disp(' '); +disp(' '); +disp('appuyez sur une touche pour continuer'); +pause; + +% calcul avec deux valeurs égales et une distincte +% premier calcul +disp(' '); +disp(' '); +disp('premier calcul avec deux valeurs égales et une autre distinctes (x2=x1 et x0~=x1)'); +disp('valeurs des trois différences divisées f[x0],f[x0,x1],f[x0,x0,x1]'); +disp(' '); +D1A=simplify(limit(D,x2,x1)); +D1B=simplify(diff_div([x0 x1 x1],[y0 y1 subs(fp,x,x1)])); +disp('calcul par passage à la limite'); +pretty(D1A); +disp('calcul par utilisation de diff_div'); +pretty(D1B); +disp('comparaison entre les deux expressions'); +disp(simplify(D1A-D1B)); +disp(' '); +disp(' '); +disp('appuyez sur une touche pour continuer'); +pause; + +% second calcul +disp(' '); +disp(' '); +disp('second calcul avec deux valeurs égales et une autre distinctes (x0=x1 et x2~=x1)'); +disp('valeurs des trois différences divisées f[x0],f[x0,x0],f[x0,x0,x2]'); +disp(' '); +D1C=simplify(limit(D,x1,x0)); +D1D=simplify(diff_div([x0 x0 x2],[y0 subs(fp,x,x0) y2])); +disp('calcul par passage à la limite'); +pretty(D1C); +disp('calcul par utilisation de diff_div'); +pretty(D1D); +disp('comparaison entre les deux expressions'); +disp(simplify(D1C-D1D)); +disp(' '); +disp(' '); +disp('appuyez sur une touche pour continuer'); +pause; + +% calcul avec les trois valeurs égales +disp(' '); +disp(' '); +disp('calcul avec les valeurs égales (x0=x1=x2)'); +disp('valeurs des trois différences divisées f[x0],f[x0,x0],f[x0,x0,x0]'); +disp(' '); +D2A=simplify(limit(limit(D,x2,x1),x1,x0)); +D2B=simplify(diff_div([x0 x0 x0],[y0 subs(fp,x,x0) subs(fpp,x,x0)])); +disp('calcul par passage à la limite'); +pretty(D2A); +disp('calcul par utilisation de diff_div'); +pretty(D2B); +disp('comparaison entre les deux expressions'); +disp(simplify(D2A-D2B)); +disp(' '); +disp(' '); +disp('appuyez sur une touche pour continuer'); +pause; + +% un calcul de dérivée +disp(' '); +disp(' '); +disp('calcul de la derivee de F : x-> f[x0,x1,x]'); +disp(' '); +dodo=diff_div([x0 x1 x],[y0 y1 f]); +F=dodo(3); +D3A=simplify(diff(F,x)); +dudu=diff_div([x0 x1 x x],[y0 y1 f fp]); +D3B=simplify(dudu(4)); +disp('expression de F''(x)'); +pretty(D3A); +disp('appuyez sur une touche pour continuer'); +pause; +disp('expression de f[x0,x1,x,x]'); +pretty(D3B); +disp(' '); +disp('comparaison entre les deux expressions'); +disp(simplify(D3A-D3B)); +disp(' '); +disp(' '); +disp('appuyez sur une touche pour continuer'); +pause; + +% dernier calcul : erreur !! +disp('essai de calcul de f[x0],f[x0,x1],f[x0,x1,x1] par substitution maladroite'); +pretty(subs(D,x2,x1)); \ No newline at end of file diff --git a/tel/TPmatlab/interpolation/TP2f/demo_trace_diff_div.m b/tel/TPmatlab/interpolation/TP2f/demo_trace_diff_div.m new file mode 100644 index 0000000..daa1919 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/demo_trace_diff_div.m @@ -0,0 +1,73 @@ +% demo_trace_diff_div : script tracé de polynômes d'interpolation pour des supports quelconques. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div, eval_horner +% +% ********************************************************* + +clear all; +choix=input('entrer les points à la main (taper m) ou en vecteur (taper v) : ','s'); +if strcmp(choix,'m') + n=input('Entrez l''entier n strictement positif (n-1 est le nombre total du point du support, en comptant les répétitions) : '); + auxix=input('Entrez le point du support numéro 0 : '); + auxiy=input('Entrez la valeur de la fonction en ce point : '); + Xb(1)=auxix; + Yb(1)=auxiy; + X=zeros(1,n+1); + Y=zeros(1,n+1); + X(1)=auxix; + Y(1)=auxiy; + r=auxix; + k=0; + for i=1:n + auxix=input(['Entrez le point du support numéro ',int2str(i),' : ']); + X(i+1)=auxix; + if (X(i+1))~=r + r=X(i+1); + k=0; + Xb=[Xb auxix]; + else + k=k+1; + end + if (k==0) + auxiy=input('Entrez la valeur de la fonction en ce point : '); + Y(i+1)=auxiy; + Yb=[Yb auxiy]; + else + Y(i+1)=input(['Entrez la valeur de la dérivée ',int2str(k),... + '-ième de la fonction en ce point : ']); + end + end +else + X=input('Entrez le vecteur des abscisses : '); + Y=input('Entrez le vecteur des ordonnées : '); + n=length(X)-1; + Xb(1)=X(1); + Yb(1)=Y(1); + r=X(1); + for i=1:n + auxix=X(i+1); + if (auxix)~=r + r=auxix; + Xb=[Xb auxix]; + Yb=[Yb Y(i+1)]; + end + end +end +pgra=input('Entrez le nombre de points pour le graphique : '); +a=X(1); +b=X(n+1); +xgra=a:(b-a)/pgra:b; +D=diff_div(X,Y); +ygra=eval_horner(xgra,X(1:end-1),D); +hold on; +plot(Xb,Yb,'go',xgra,ygra); +legend('Points du support','polynôme interpolateur'); +title(['interpolation avec ',int2str(n+1),' points']); +hold off; + + + + diff --git a/tel/TPmatlab/interpolation/TP2f/demo_verifie_formule.m b/tel/TPmatlab/interpolation/TP2f/demo_verifie_formule.m new file mode 100644 index 0000000..3df0bdf --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/demo_verifie_formule.m @@ -0,0 +1,60 @@ +% script demo_verifie_formule : verification formelle de différences divisées. +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div, eval_horner +% +% ********************************************************* + + +clear all; +syms x a b fa fpa fb fpb; +X=[a a b b]; +Y=[fa fpa fb fpb]; +D=diff_div(X,Y); +disp('affichage des différentes différences divisées : '); +disp(' '); +disp('Ici, fpa et fpb désignent f''(a) et f''(b)'); +disp(' '); +disp('f[a] : '); +disp(D(1)); +disp('f[a,a] : '); +disp(D(2)); +disp('f[a,a,b] : '); +disp(D(3)); +disp('f[a,a,b,b] : '); +disp(D(4)); +disp('appuyez sur une touche pour continuer '); +pause; +disp(' '); +disp(' '); +disp('vérification des formules : '); +disp('f[a]=f(a) : '); +auxi=simplify(D(1)-(fa)); +if (auxi==0) + disp('ok !!!'); +else + disp('erreur !!!'); +end +disp('f[a,a]=f''(a) : '); +auxi=simplify(D(2)-(fpa)); +if (auxi==0) + disp('ok !!!'); +else + disp('erreur !!!'); +end +disp('f[a,a,b]=(f(b)-f(a))/((b-a)^2)-(f''(a))/(b-a) : '); +auxi=simplify(D(3)-((fb-fa)/((b-a)^2)-(fpa)/(b-a))); +if (auxi==0) + disp('ok !!!'); +else + disp('erreur !!!'); +end +disp('f[a,a,b,b]=(f''(a)+f''(b))/((b-a)^2)+2(f(a)-f(b))/((b-a)^3) : '); +auxi=simplify(D(4)-((fpa+fpb)/((b-a)^2)+2*(fa-fb)/((b-a)^3))); +if (auxi==0) + disp('ok !!!'); +else + disp('erreur !!!'); +end diff --git a/tel/TPmatlab/interpolation/TP2f/diff_div.m b/tel/TPmatlab/interpolation/TP2f/diff_div.m new file mode 100644 index 0000000..65c0ca4 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/diff_div.m @@ -0,0 +1,70 @@ +function D=diff_div(X,Y) + + +% diff_div : calcul des différences divisées pour des points de supports quelconques +% +% ********************************************************* +% +% D=diff_div(X,Y) +% les différences divisées sont calculées par l'algorithme +% pyramidal donné en TD. +% +% variables d'entrées : +% * X : contient les valeurs x_i, pour 0 <=i<=n (dans l'ordre croissant) +% * Y : contient les valeurs f^(r)(x_i), pour 0 <=i<=n, de telle sorte que : +% si x_k est multiple avec x_{k-1}~=x_k=x_{k+1}=...=x_{k+r}~=x_{k+r+1} alors +% y_k=f(x_k), y_{k+1}=f'(x_k), ..., y_{k+r}=f^(r)(x_k). +% +% variables de sortie +% * D : contient les différences divisées (généralisée à support quelconque) +% f[x_0], f[x_0,x_1], ...., f[x_0,x_1,...,x_n] +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% val_derivee_multiple +% +% ********************************************************* + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=2 + error('nombre d''arguments de la fonction incorrect'); +end +% taille des deux variables +n=length(X)-1; +np=length(Y)-1; +if n~=np + error('les deux tableaux n''ont pas la même taille'); +end +% vérification de l'ordre croissant des x_i +% (si X non formel) +if isnumeric(X) + if max(abs(sort(X)-X))~=0 + error('support non ordonné'); + end +end + +% Corps d'algorithme +i=0; +while (i<=n) + [alpha,valfder]=val_derivee_multiple(i,X,Y); + D(i+1:i+alpha)=valfder(1)*ones(1,alpha); + i=i+alpha; +end +factoi=1; +for i=1:n; + factoi=factoi*i; + for j=n:-1:i; + dx=X(j+1)-X(j-i+1); + if (dx==0) + [alpha,valfder]=val_derivee_multiple(j-i,X,Y); + D(j+1)=(valfder(i+1))/factoi; + else + D(j+1)=(D(j+1)-D(j))/(dx); + end + end +end diff --git a/tel/TPmatlab/interpolation/TP2f/test_diff_div.m b/tel/TPmatlab/interpolation/TP2f/test_diff_div.m new file mode 100644 index 0000000..6a43a50 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/test_diff_div.m @@ -0,0 +1,34 @@ +function R=test_diff_div(T,X,Y) + +% test_diff_div : fonction test pour le calcul du polynôme interpolateur support quelconque +% +% +% ********************************************************* +% +% R=test_diff_div(T,X,Y) +% calcul de la valeur du polynôme d'interpolation défini par le nuage de point X,Y +% sur le tableau T (utilisation du calcul des différences divisées et +% de l'algorithme d'évaluation de Horner) pour support quelconque +% +% variables d'entrées : +% * T : contient le tableau rectangulaire des réels où est évalué p_n +% * X : contient les valeurs x_i, pour 0 <=i<=n (dans l'ordre croissant) +% * Y : contient les valeurs f^(r)(x_i), pour 0 <=i<=n, de telle sorte que : +% si x_k est multiple avec x_{k-1}~=x_k=x_{k+1}=...=x_{k+r}~=x_{k+r+1} alors +% y_k=f(x_k), y_{k+1}=f'(x_k), ..., y_{k+r}=f^(r)(x_k). +% +% variables de sortie +% * R : contient le tableau rectangulaire des images de T par p_n, +% polynôme d'interpolation de f sur le support x_i, pour 0 <=i<=n +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div, eval_horner +% +% ********************************************************* + + +D=diff_div(X,Y); +R=eval_horner(T,X(1:end-1),D); + diff --git a/tel/TPmatlab/interpolation/TP2f/val_derivee_multiple.m b/tel/TPmatlab/interpolation/TP2f/val_derivee_multiple.m new file mode 100644 index 0000000..f67665b --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/val_derivee_multiple.m @@ -0,0 +1,64 @@ +function [alpha,valfder]=val_derivee_multiple(i,X,Y) + + +% val_derivee_multiple : sortie des valeurs des dérivées successives pour un support quelconque +% +% ********************************************************* +% +% [alpha,valfder]=val_derivee_multiple(i,X,Y) +% sorties des valeurs des dérivées successives (définies dans Y) +% pour un support à n+1 points quelconque (ordonnée) stocké dans X, +% pour un élément x_i de ce support. +% +% +% variables d'entrées : +% * i : entier compris entre 0 et n +% * X : contient les valeurs x_j, pour 0 <=j<=n (dans l'ordre croissant) +% * Y : contient les valeurs f^(r)(x_j), pour 0 <=j<=n, de telle sorte que : +% si x_k est multiple avec x_{k-1}~=x_k=x_{k+1}=...=x_{k+r}~=x_{k+r+1} alors +% y_k=f(x_k), y_{k+1}=f'(x_k), ..., y_{k+r}=f^(r)(x_k). +% +% variables de sortie +% * alpha : nombre d'éléments du n+1 uplet (x_0,x_1,...,x_n) égaux à x_i. +% * valfder : vecteur ligne de longueur alpha tel que : +% pour tout l dans {0,...,alpha-1}; valfder(l+1)=f^{l}(x_i). +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% aucune +% +% ********************************************************* +% + + +% Contrôles d'entrée +% +% nombre d'arguments +if nargin~=3 + error('nombre d''arguments de la fonction incorrect'); +end +% taille des deux variables +n=length(X)-1; +np=length(Y)-1; +if n~=np + error('les deux tableaux n''ont pas la même taille'); +end +% vérification de l'ordre croissant des x_i +% (si X non formel) +if isnumeric(X) + if max(abs(sort(X)-X))~=0 + error('support non ordonné'); + end +end +% vérification de l'appartenance de i à {0,...,n} +% vérification de l'ordre croissant des x_i +if ~((0<=i) & (i<=n)) + error('le premier argument n''appartient pas à {0,...,n}'); +end + + +% Corps d'algorithme +k=find(X==X(i+1)); +alpha=length(k); +valfder=Y(k); \ No newline at end of file diff --git a/tel/TPmatlab/interpolation/TP2f/verifie_diff_div.m b/tel/TPmatlab/interpolation/TP2f/verifie_diff_div.m new file mode 100644 index 0000000..902c7b4 --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/verifie_diff_div.m @@ -0,0 +1,53 @@ +function res=verifie_diff_div(X,Y) + +% verifie_diff_div : verification pour le calcul du polynôme interpol. support quelconque +% +% +% ********************************************************* +% +% res=verifie_diff_div(X,Y) +% vérification du calcul du polynôme d'interpolation défini par le nuage de point X,Y +% pour support quelconque +% +% variables d'entrées : +% * X : contient les valeurs x_i, pour 0 <=i<=n (dans l'ordre croissant) +% * Y : contient les valeurs f^(r)(x_i), pour 0 <=i<=n, de telle sorte que : +% si x_k est multiple avec x_{k-1}~=x_k=x_{k+1}=...=x_{k+r}~=x_{k+r+1} alors +% y_k=f(x_k), y_{k+1}=f'(x_k), ..., y_{k+r}=f^(r)(x_k). +% +% variables de sortie +% * res : contient le maximum des valeur absolue de +% n-p où p est le degré du polynôme d'interpolation p_n +% et n+1 est le nombre de points. +% maximum sur i dans {0,...,n} maximum sur l dans alpha_i-1 +% de la quantitée p_n^(l)(x_i)-f^(l)(x_i) (dérivée l fois) où alpha_i est +% le nombre d'éléments du n=1-uplet(x_0,...,x_n) égaux à x_i. +% +% Le calcul est correct si et seulement si res est nul. +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div, newton_cano, val_derivee_multiple +% +% ********************************************************* + +D=diff_div(X,Y); +P=newton_cano(D,X(1:end-1)); +n=length(X)-1; +p=length(P)-1; +erdeg=abs(n-p); +er=0; +i=0; +while (i<=n) + [alpha,valfder]=val_derivee_multiple(i,X,Y); + Q=P; + er=max(er,abs(polyval(P,X(i+1))-valfder(1))); + for l=1:alpha-1 + Q=polyder(Q); + er=max(er,abs(polyval(Q,X(i+1))-valfder(l+1))); + end + i=i+alpha; +end +res=max(er,erdeg); diff --git a/tel/TPmatlab/interpolation/TP2f/verifie_diff_div_symb.m b/tel/TPmatlab/interpolation/TP2f/verifie_diff_div_symb.m new file mode 100644 index 0000000..a5252ac --- /dev/null +++ b/tel/TPmatlab/interpolation/TP2f/verifie_diff_div_symb.m @@ -0,0 +1,57 @@ +function res=verifie_diff_div_symb(X,Y) + +% verifie_diff_div_symb : verification pour le calcul du polynôme interpol. support quelconque en symbolique +% +% +% ********************************************************* +% +% res=verifie_diff_div_symb(X,Y) +% vérification du calcul du polynôme d'interpolation défini par le nuage de point X,Y +% pour support quelconque où X et Y sont des vecteurs symboliques. +% Vérification en formel. +% +% +% variables d'entrées : +% * X : contient les valeurs x_i, pour 0 <=i<=n (dans l'ordre croissant) +% * Y : contient les valeurs f^(r)(x_i), pour 0 <=i<=n, de telle sorte que : +% si x_k est multiple avec x_{k-1}~=x_k=x_{k+1}=...=x_{k+r}~=x_{k+r+1} alors +% y_k=f(x_k), y_{k+1}=f'(x_k), ..., y_{k+r}=f^(r)(x_k). +% +% variables de sortie +% * res : contient le maximum des valeur absolue de +% n-p où p est le degré du polynôme d'interpolation p_n +% et n+1 est le nombre de points. +% maximum sur i dans {0,...,n} maximum sur l dans alpha_i-1 +% de la quantitée p_n^(l)(x_i)-f^(l)(x_i) (dérivée l fois) où alpha_i est +% le nombre d'éléments du n=1-uplet(x_0,...,x_n) égaux à x_i. +% +% Le calcul est correct si et seulement si res est nul +% et incorrect si res est egal à un +% +% +% +% ************ Fonctions auxiliaires utilisées ************ +% +% diff_div, eval_horner, val_derivee_multiple +% +% ********************************************************* + +D=diff_div(X,Y); +syms x; +R=eval_horner(x,X(1:end-1),D); +er=0; +i=0; +n=length(X)-1; +while (i<=n) + [alpha,valfder]=val_derivee_multiple(i,X,Y); + Q=R; + erloc=~(simplify(subs(R,'x',X(i+1)))==valfder(1)); + er=max(er,erloc); + for l=1:alpha-1 + Q=diff(Q); + erloc=~(simplify(subs(Q,'x',X(i+1)))==valfder(l+1)); + er=max(er,erloc); + end + i=i+alpha; +end +res=er; \ No newline at end of file