Commit 9bc6e052 authored by Conor McCoid's avatar Conor McCoid
Browse files

Assistant work: added TP5 from 2020

parent 390a39de
function [y,fe]=Adapt(f,a,b,abstol,fe)
% ADAPT should not be called by itself, it is used by Integrate
% [y,fe]=Adapt(f,a,b,abstol) integrates f in the interval [a,b].
% It is the estimate of the integral using Monte Carlo, abstol is the
% absolute tolerance, and fe is the number of function evaluations needed.
I1=Simpson(f,a,b);
I2=Gauss(f,a,b,5);
err=abs(I1-I2);
if(err<=abstol)
y=I2; fe=fe+5;
else
m=(a+b)/2;
[y1,fe]=Adapt(f,a,m,abstol,fe);
[y2,fe]=Adapt(f,m,b,abstol,fe);
y=y1+y2;
end
\ No newline at end of file
function [I]=Gauss(f,a,b,n)
[pts,w]=GaussCoefficients(n);
h=b-a;
I=h*sum(w.*f(a+h*pts));
\ No newline at end of file
function [pts,w]=GaussCoefficients(s)
for i=1:s-1
gamma(i)=i/sqrt((2*i-1)*(2*i+1));
end
A=diag(gamma,-1)+diag(gamma,1);
[W,P]=eig(A);
pts=diag(P); pts=(1+pts')/2;
w=W(1,:).^2;
\ No newline at end of file
function [y,fe] = Integrate(f,a,b,tol)
%INTEGRATE Adaptive integration of f using Simpson's rule and Gauss' rule with 5 nodes
% [y,fe]=Integrate(f,a,b,tol) integrates f in the interval [a,b] using
% Simpson and Gauss' rule. The relative tolerance is specified by tol. fe returns
% the number of function evaluations needed.
%Verify tol
if(tol>=1)
disp('Choose a smaller tolerance');
elseif(tol<eps)
disp('Choose a bigger tolerance');
else
%Monte Carlo estimation
x=rand(1,5)*(b-a) +a; m=(a+b)/2;
I=(b-a)/8*(f(a)+f(b)+f(m)+sum(f(x)));
[y,fe]=Adapt(f,a,b,abs(I)*tol,0);
end
function [I]=Simpson(f,a,b)
m=(a+b)/2;
I=(b-a)/6*(f(a)+4*f(m)+f(b));
\ No newline at end of file
f1=@(x) 4./(1+x.^2);
f2=@(x) 4.*sqrt(1-x.^2);
for i=1:15
tol(i)=10^(-i);
[y1(i),fe1(i)] = Integrate(f1,0,1,tol(i));
err1(i)=abs(pi-y1(i));
[y2(i),fe2(i)] = Integrate(f2,0,1,tol(i));
err2(i)=abs(pi-y2(i));
end
figure(1)
semilogy(fe1,tol,'b-');
hold on
semilogy(fe1,err1,'r.-');
title('f=4/(1+x^2)')
legend('Estimated error','Real error')
figure(2)
semilogy(fe2,tol,'b.-');
hold on
semilogy(fe2,err2,'r.-');
title('f=4*sqrt(1-x^2)');
legend('Estimated error','Real error')
\ No newline at end of file
\documentclass[12pt,a4paper]{article}
\usepackage{amsmath} % Pour les symboles complementaire comme les matrices !
\usepackage{amssymb,amsfonts}
\usepackage{alltt}
\usepackage{epsfig}
\usepackage[french]{babel}
%\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
%\usepackage[T1]{fontenc} % Pour la bonne cesure du francais
%\usepackage[applemac]{inputenc}
%\usepackage{inputenc}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{example}{Example}
\newcommand{\qed}{\hfill$\qedsquare$\goodbreak\bigskip}
\advance\voffset by -30mm \advance\hoffset by -25mm
\setlength{\textwidth}{175mm} \setlength{\textheight}{260mm}
\pagestyle{empty}
\newdimen\figcenter
\newdimen\figlarge
\newdimen\figheight
\def\figinsert#1{
\figcenter=\hsize\relax
\advance\figcenter by -\figlarge
\divide\figcenter by 4
\vskip +1.truecm
\vskip\figheight\noindent\hskip\figcenter
\special{psfile=#1}
\vskip -.7truecm}
\renewcommand{\vec}[1]{\mathbf{#1}}
\begin{document}
\begin{center}\noindent {\large UNIVERSIT\'E DE GEN\`EVE} \end{center}
\noindent Section de Math\'ematiques \hfill\`A rendre le jeudi 28 novembre 2019.\\
\noindent Facult\'e des Sciences \hfill L'évaluation aura lieu le vendredi 29 novembre 2019.\\[1mm]
\hrule
\bigskip
\bigskip
\begin{center} {\bf {\Large Analyse Num\'erique}}\end{center}
\begin{center} {\bf {\Large Travaux Pratiques~-~S\'erie 5}}\end{center}
\bigskip
\begin{enumerate}
\item La règle du trapèze composée donne des résultats étonnamment bons quand on l'applique aux fonctions périodiques (i.e., $f^{(k)}(a) = f^{(k)}(b)$, pour
$k=0,1,\dots$.), bien meilleures que l'estimation d'erreur vue en cours pour des fonctions générales.
\begin{enumerate}
\item Intégrer la fonction périodique
\[
\mathcal I = \int_0^1\exp\big(x^2(1-x)^2\big) dx
\]
en utilisant la règle du trapèze composée
(avec des n\oe{}uds équidistants $a=x_0 < x_1<\dots<x_n= b$).
Tracer l'erreur d'intégration comme fonction du nombre de n\oe{}uds $n=2:2:52$ (voir Fig.~\ref{fig:period}), où vous pouvez utiliser la fonction de MATLAB \texttt{integral} pour calculer $I$ `exactement'. Comparer les résultats avec les intégrales obtenues par la méthode de Simpson composée.
Répéter l'expérience pour la fonction non périodique $\int_0^1 \exp(-x^2) dx$ pour voir la différence.
\item Expliquer (non rigoureusement) pourquoi la règle du trapèze fonctionne aussi bien.
{\bf Indice:}
(1) La règle du trapèze est exacte pour les fonctions trigonométriques $\sin(2k\pi x)$ et $\cos(2k\pi x)$ sur $x\in[0,1]$, pour $k\in \mathbb Z$. Vérifier cela numériquement pour quelques valeurs $k$.
(2) Les fonctions périodiques admettent une extension en séries de Fourier.
\end{enumerate}
\item Le but de cet exercice est d'implémenter un intégrateur \emph{adaptatif}
pour évaluer des intégrales du type
$
\int_a^b f(x)\,dx.
$
Le principe d'un programme adaptatif est de raffiner le découpage de l'intervalle si l'erreur
locale n'est pas suffisamment petite.
\begin{enumerate}
\item \'Ecrire une programme pour calculer les poids et les n\oe uds de la formule de Gauss à $s$~n\oe uds
en diagonalisant une certaine matrice tridiagonale (voir le résultat de l'exercice 2, série 9). Utiliser l'en-tête suivant :
\small
\begin{verbatim}
function [c,b] = GaussCoefficients(s)
% GAUSSCOEFFICIENTS Computations of nodes and weights for s-nodes Gauss' rule
% [c,b] = GaussCoefficients(s) computes the nodes and the weights for
% the s-nodes Gauss' rule in the interval [0,1].
\end{verbatim}
Tester votre programme en cherchant les n\oe uds et les poids pour $s=4$ et $s=5$. \\
Vous devriez obtenir pour $s=4$ :
$$c_1\approx 0.0694\, , \,c_2\approx 0.33 \,,\, c_3\approx0.67 \,,\, c_4\approx 0.9306\; ;$$$$ \; b_1 = b_4 \approx 0.1739 \,,\, b_2 = b_3 \approx 0.3261.$$
et pour $s=5$ :
$$c_1\approx 0.0469\, , \,c_2 \approx 0.2308 \,,\, c_3\approx 0.5 \,,\, c_4 \approx 0.7692\,,\, c_5\approx 0.9531\; ;$$$$ \; b_1 = b_5 \approx 0.1185\,,\, b_2 = b_4 \approx 0.2393 \,,\, b_3 \approx 0.2844.$$
\item Dans cette deuxième partie, on va écrire un programme adaptatif
pour évaluer les intégrales numériquement. Les formules de quadrature que nous allons utiliser sont la règle de Simpson:
$$ \int_a^b f(x)dx \approx \frac{b-a}{6}(f(a) + 4f(m) + f(b)),\;\;\;\;\;\mbox{}\; m = \frac{1}{2}(a+b),$$
et la règle de Gauss à 5 n\oe uds :
$$\int_a^b f(x)\,dx = h\,\displaystyle\sum_{i=1}^s\,b_if(a+hc_i),\;\;\;\;\;\mbox{}\; h=(b-a).$$
(utiliser les coefficients trouvés dans l'exercice 1 avec autant de précision que possible (16 chiffres décimaux)).
\begin{enumerate}
\item (4) Commencer votre programme avec l'en-tête suivant:
\small
\begin{verbatim}
function [y,fe] = Integrate(f,a,b,tol)
%INTEGRATE Adaptive integration of f using Simpson's rule and Gauss' rule with 5 nodes
% [y,fe]=Integrate(f,a,b,tol) integrates f in the interval [a,b] using
% Simpson and Gauss' rule. The relative tolerance is specified by tol. fe returns
% the number of function evaluations needed.
\end{verbatim}
\normalsize
Votre fonction {\tt Integrate} doit exécuter les tâches suivantes:
\begin{itemize}
\item Vérifier que la tolérance {\tt tol} est raisonnable
(p. ex. entre {\tt eps} et $10^{-4}$);
\item Obtenir une première estimation $I_s$ de l'intégrale à l'aide
d'une méthode de Monte Carlo à cinq points:
$$ I_s := \frac{b-a}{8}\left(f(a)+f(b)+f(m)+\sum_{j=1}^5 f(x_j)\right),$$
$x_1,\dots,x_5$ sont cinq points aléatoires dans l'intervalle $[a,b]$ et $m = \frac{1}{2}(a+b)$ le milieu de l'intervalle;
\item Appeler la fonction interne {\tt Adapt}, que vous écrirez \`a
l'étape suivante.
\end{itemize}
\item (6) La fonction interne {\tt Adapt} est définie de la manière suivante:
\small
\begin{verbatim}
function [y,fe]=Adapt(f,a,b,abstol)
% ADAPT should not be called by itself, it is used by Integrate
% [y,fe]=Adapt(f,a,b,abstol) integrates f in the interval [a,b].
% It is the estimate of the integral using Monte Carlo, abstol is the
% absolute tolerance, and fe is the number of function evaluations needed.
\end{verbatim}
\normalsize
Pour estimer l'erreur locale, on considère la différence
$err = |I_1 - I_2|$ des approximations ci-dessous, et on compare cette différence
à la tolérance absolue $abstol = |I_s| * tol$.
\begin{itemize}
\item $I_1$: Le résultat de la règle de Simpson sur l'intervalle $[a,b]$,
\item $I_2$: Le résultat de la règle de Gauss à 5 n\oe uds sur l'intervalle $[a,b]$.
\end{itemize}
Si $err \leq abstol$, la fonction retourne $y :=I_2 $.
Si $err > abstol$, la fonction subdivise l'intervalle $[a,b]$ en $[a,m]\cup[m,b]$ et s'appelle elle-même
sur chacun des sous-intervalles.
\item (4) Faire un graphe en échelle logarithmique de l'erreur en fonction du nombre
d'évaluations $fe$ pour les deux estimations de $\pi$ suivantes:
$$\pi = \int_0^1 \frac{4\,dx}{1+x^2},\qquad \pi = \int_0^14\sqrt{1-x^2}dx.$$
\end{enumerate}
\end{enumerate}
\end{enumerate}
\begin{figure}[h]
\begin{center}
\includegraphics[width=.45\textwidth]{p1.pdf}~\includegraphics[width=.45\textwidth]{p2.pdf}
\end{center}
\caption{Graphique de l'exercice 1:
(Gauche)
$f(x) = \exp(x^2(1-x)^2)$;
(Droite)
$f(x) = \exp(-x^2)$.
}
\label{fig:period}
\end{figure}
%\begin{center}
%---------------------------------
%\end{center}
\end{document}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment