TRAVAUX PARATIQUES

D’ANALYSE NUMERIQUE
 Initiation à la programmation
 Travaux pratiques d’analyse numérique
 Petits projets d’application

Abdelaziz OUATIZERGA
Docteur en physique à l’université Houari Boumediene – Bab Ezzouar

Djafer Benrabia
Docteur en physique à l’université Houari Boumediene – Bab Ezzouar
USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Partie I :
Initiation à la programmation
 Prise en main d’Octave
 Travailler avec la fenetre de commandes
 Ecriture de scripts simples
 Ecriture de scripts évolués
o L’instruction IF
o Les boucles For, while
o Les fonctions
 Ecriture des scripts de fonctions

1
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Prise en main du logicel Octave

1 : Le menu

2 : La barre d’outils

3 : La fénêtre des commandes

4 : L’éditeur de scripts

La fenêtre de commande (command window) est la fentre la plus utilisée. Elle permet d’exécuter les
commandes de tracé de courbe, d’éxécuter des commandes de calcul entre variable d’éxécuter des fonctions
intégré, des fonctions personnelles crées par l’utilisateur et des scripts de commandes .etc

Toute commande dans cette fenetre est exécuter en appyant sur la touche entrée.

Exemple :

>> x=1.54

2
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Après avoir tapé sur entrée, la réponse suivante est affiché.

x = 1.5400
La commande précédente affecte la valeur 1.54 à la variable x.

A chaque exécution d’une commande, il y a une réponse de la sorte. Lorsqu’on a afaire à plusieurs
commandes et que l’affichage de ceraines commandes dvient enuyeux ; on peut désactiver l’affichage en
terminant ces commandes par un point virgule.

La fnétre de commande est accessible à partir du menu principal : Fenetres / Fenetre des commandes

4 : L’éditeur de script

3
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Initiation à la programmation
ave Octave
Travailler avec la fenetre de commandes
Dans ce chapitre, nous allons détailler les principales fonctions d’Octave utilisées dans les séances des tra-
vaux pratiques. Chacune des commandes en gris sont à retaper sur la fenêtre de commande d’Octave de votre
ordinateur.
Certaines remarques sont à prendre en compte :

 Les résultats des commandes ne sont pas affichés intentionnellement. En effet, c’est un document
voulu être pratique qui ne sera efficace qu’avec l’effort de l’étudiant.
 Le caractère % (pourcentage) sert à commenter le texte qui le suit, donc n’est pas pris en compte par
Octave et par conséquent n’est pas à retaper obligatoirement.
 Les commandes fonctionnent indifféremment sur le logiciel libre « Octave » et sur le logiciel com-
mercial « Matlab ».

Pour exécuter MATLAB / OCTAVE il faut faire un clic sur son icône dans le menu principal du système..

Vous pouvez quitter par la ligne de commande en tapant : quit

Recherche de l’aide :

Lors du travail sur Octave, nous sommes amenés maintes fois à chercher une description dé-
taillée d’une commande ou une d’une certaine fonction ou un mot clé. Pour ce faire, les
commandes suivantes peuvent être utilisées sur la fenêtre commandes.

≫ help % ou à partir à partir du menu

≫ lookfor % dans la fenêtre de commande
≫ demo matlab % pour apprendre facilement (Matlab uniquement)

Les opérations du calcul

Toutes les opérations disponibles dans une calculatrice sont réalisables dans la fenêtre de com-
mandes d’Octave avec les opérateurs du calcul habituels tels que : + - * / ^
La virgule décimale est représentée exclusivement par un point s Octave. La virgule étant utilisée dans Oc-
tave pour séparer les paramètres et les commandes.

>> 1+2 % addition

>> 5*3 % multiplication
>> 2.5^5 % puissance → (2.5)5
>> sqrt (3.256) % racine carré → √3,256
>> log (30) % logarithme népérien → ln(30)
>> log10 (30) % logarithme base 10 → log(30)

4
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

>> cos(1.2) % cosinus de l’angle 1.2 radian

>> sind(30) % sinus de 30 degré
>> exp(7.5) % exponentielle → 𝑒 7.5
>> 1e-10 % 1e3 =103 , 51e-3 = 0.51, 3e8=3x108 ... etc.

La commande format permet de modifier la précision de l’affichage du résultat :

Remarquer le format du résultat :

≫ format long
≫ 7 + 10/3 + 5 ^ 1.2
≫ format short
≫ 7 + 10/3 + 5 ^ 1.2

Les variables
Les variables servent à stocker des valeurs numériques ou des chaines de caractères pour une utilisation ulté-
rieure. Les noms de variables doivent respecter les critères suivants :
 Ne doivent pas contenir des espaces
 Ne doivent pas contenir des caractères spéciaux tel que : / ) + - * % :
 Ne doivent pas commencer par un chiffre.

Important: Octave respecte la casse c’est-à-dire une lettre majuscule est différente de la lettre
minuscule.
La variable x , par exemple, est différente de la variable X.

Dans Octave, si aucune variable n’est spécifiée, la réponse de chaque opération de calcul est
stockée par défaut dans la variable ans
La variable ans peut être utilisé par la suite dans d’autres opérations.

≫2+3
≫x= ans +5/ans

On peut spécifier librement le nom d’une variable à condition de respecter les conditions précé-
dentes.

Taper les commandes suivantes et voir le résultat à chaque fois. Utilisez des noms de variables de
votre choix.

≫x=3 ; % affecter la valeur 3 à la variable nommée x

≫X = x+5 ; % Calculer x+3 puis affecter le résultat à X.
≫z = sqrt(X) ; % Calculer la racine carré de X et affecter à z
≫disp(z) % Afficher la valeur de z

≫abc="salam alaikom" ; % abc est de type chaine de caractère

≫abc % afficher le contenu de la variable abc

Les variables sont de type entiers, réels, chaine de caractères ou logiques

Dans Octave, le type d’une variable est déterminé par la valeur affectée à cette variable

5
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Donc le type de la variable change si on change son contenu

≫ abc = 1250 % la variable abc devient de type entier

Pour voir les attributs des variables présents dans l’espace de travail d’Octave :

≫ who % afficher les variables dans la mémoire

≫ whos % afficher plus de détails sur ces variables

≫ clear x % effacer la variable x

≫ clear % effacer toutes les variables de la mémoire
≫ clear all % effacer toutes les variables de la mémoire

Remarque :

Les opérations peuvent être écrites sur une même ligne de commande, séparées par des points vir-
gules.
Le point-virgule sert à délimiter une commande. Lorsqu’elle est présente, le résultat de la com-
mande ne s’affiche pas.

≫var1= 2 ; tmp = var1 + 2.5 ; z = tmp^2.3 ;

≫z % afficher la valeur de z

Exercices
Exercice 1 : Trouver les deux solutions de l’équation au second degré suivante :

183.25𝑥 2 + 17.33𝑥 − 389.15 = 0

Suivre les étapes suivantes ::

1- Affecter les coefficients de l’équation à des variables de votre choix. Soit par exemple a, b, c ces va-
riables
2- Ecrire l’expression qui permet de calculer le discriminant. Soit delta le nom de la variable
3- Ecrire les expressions des solutions. Soit x1 et x2 les noms des deux variables
4- Afficher les résultats

Exercice 2 : Sans retaper toutes les expressions résoudre l’équation :

13.5𝑥 2 − 30.05𝑥 − 81.2 = 0

Conseil : Utiliser la flèche du haut du clavier pour rappeler les expressions de a, b, c, delta, x1 et x2.

Les vecteurs
Un vecteur est un tableau à une dimension. C’est une liste de valeurs de même type désignée par un
nom de variable qu’on définit au choix.
Dans Octave, cette liste de valeurs est appelée vecteur. Les valeurs sont appelés les éléments du
vecteur.

6
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Les éléments du vecteur sont listé entre crochets [ . . ]

Séparateurs : blanc, virgule ou point-virgule.

≫ v= [1 3 1 -6 8]

Le vecteur v a quatre éléments de type entier.

Vecteurs lignes : le séparateur est une virgule ou un espace

≫ v= [1 3 1 -6 8]
≫ w= [2, -2, 1, 6, - 8]

Vecteurs colonnes : le séparateur est un point-virgule

≫ u = [7; 8; 9; 10; 11] ;

Aussi, on peut transformer un vecteur ligne en un vecteur colonne et vice-versa en prenant le transposé.

≫ u’ % afficher le transposé de u
≫ u=u’ % le vecteur u devient un vecteur colonne.
≫ u=u’ % le vecteur u devient de nouveau un vecteur ligne.

Voici quelques fonctions pratiques intégrées dans Octave sur les vecteurs :

length(u) : donne le nombre d’élément du le vecteur u

sum(u) : calcule la somme de tous les éléments du vecteur u.
prod(u) : calcule le produit de tous les éléments du vecteur u.
norm(u) : calcule la norme du vecteur u

Attention :

Avant de faire des opérations entre vecteurs, il est important de se rassurer si les deux vecteurs sont du même
type ou non. Car, si les deux vecteurs sont du même type le résultat est un vecteur sinon le résultat sera une
matrice.

Les éléments du vecteur sont numérotés à partir de 1. Le ieme élément du vecteur u est noté u(i)
L’entier i est l’index de l’élément (le numéro d’ordre de l’élément dans le vecteur).

≫u(1) % afficher le premier élément du vecteur u

≫u(3) % afficher le 3e élément du vecteur u
≫u(2:5) % afficher les éléments de 2 à 5 de u
≫length(u) % la longueur du vecteur u (= nombre d’éléments)

Les opérations d’addition et soustraction entre vecteurs sont réalisées à l’aide des opérateurs + et -
ordinaires.

≫ u = v + w % somme et différence entre deux vecteurs

≫ u = v - w % somme et différence entre deux vecteurs

Chaque élément du vecteur u est la somme de deux éléments de v et w du même index (du même
ordre)

7
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

De manière générale, dans Octave, les opérations entre vecteurs et matrices obéissent à des règles
du calcul matriciel mathématique. Par conséquent, le produit de deux vecteurs du même type n’est
pas autorisé.

≫ v * v % ERROR pour le produit matriciel

≫ u * v’ % produit scalaire : le résultat est un scalaire= u.*v
≫ u’*v % résultat est une matrice

Octave offre la possibilité d’effectuer les opérations de multiplication et de division de deux vec-
teurs de même type mais élément par élément à la manière de l’addition et soustraction. Pour ce
faire, il faut utiliser les mêmes opérateurs précédés par un point : .* pour la multiplication et ./ pour
la division, de même pour le calcul de la puissance : .^

Tester et examiner le résultat des opérations suivantes :

≫ u = v.* w % multiplication élément par élément

≫ u = v./ w % division élément par élément
≫ u = v.^3 % opération sur chaque élément
≫ u = cos(v) % opération sur chaque élément 𝑢(𝑖) = cos(𝑣(𝑖)) , 𝑖 = 1. . .4
≫ F = u.^2 +1 % F un vecteur tel que F(i)=u(i)^2+1, i=1 …4

Lorsqu’une opération est appliquée sur les éléments d’un vecteur ; le résultat sera donc un vecteur dont les
éléments seront les images des éléments du vecteur initial par l’opération appliquée.

Octave offre deux autres manières très pratiques pour créer des vecteurs à pas (intervalle) régulier :

La première : en spécifiant un pas d’incrémentation régulier :

v = [debut:pas:fin] ou sans les crochets : v = debut:pas:fin

Début: la valeur du premier élément
fin : la valeur du dernier élément
pas : le pas d’incrémentation (1 si n’est pas spécifié)

≫ v=[1:2:33] % définir un vecteur de 1 à 33 par pas de 2

Ou simplement :

≫ v=1:2:33 % définir un vecteur de 1 à 33 par pas de 2

La deuxième : Lorsque’ on préfère spécifier un nombre de points au lieu du pas d’incrémentation.

v = linspace (debut, fin, n)

Le paramètre n est le nombre de points à prendre entre les limites de l’intervalle debut et fin et qui sera égal
aussi au nombre d’éléments du vecteur v

≫ v = linspace (1, 33, 16)

8
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Les matrices
Définition d’une matrice :

On peut définir une matrice directement en spécifiant ses éléments ; espace ou virgule pour séparer les élé-
ments d’une ligne et point-virgule pour séparer les lignes.

On peut atteindre un élément spécifique :

≫ A(2,3) % c’est l’élément de la 2eme ligne, 3eme colonne ici 6

≫ A(2:3,1:2) % les éléments des lignes 2 et 3 et colonnes 1 et 2
≫ A(1:3, 2) % les lignes de 1 à 3 , colonne 2
≫ A(:,2) % Afficher la 2eme colonne
≫ A(3,:) % Afficher la 3eme ligne

Les opérations entre matrices :

Produit matriciel

Soit le vecteur:
≫ b = [1; 1; 1] ; % definir un vecteur colonne

≫ A*b % produit matriciel (matrice par vecteur colonne)

Division matricielle à gauche

≫ A\b % solution de A x = b ou b est un vecteur colonne (x=A\b)

Opérations spécial Matlab

Les opérateurs : .* ./ .^ agissent sur les éléments de la matrice comme pour les vecteurs

≫ A.*A % multiplication élément par élément

≫ A./A % division élément par élément

Certaines fonctions usuelles sur les matrices :

≫ zeros(n,m) % crée une matrice n lignes et m colonnes remplie de 0

≫ ones(n,m) % crée une matrice n lignes et m colonnes remplie de 1
≫ size(A) % la taille de la matrice [nb lignes, nb colonnes]
≫ diag(A) % le vecteur contenant les éléments de la diagonale de A

Réalisation des graphes

9
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Pour tracer le graphe d’une fonction 𝑓(𝑥) sur l’intervalle [𝑎, 𝑏] , on doit créer deux vecteurs nommés X et Y
par exemple. X un vecteur contenant les valeurs de l’abscisse x allant de a à b par pas régulier (mais ce n’est
pas obligatoire). Y un vecteur dont les éléments sont les images des éléments du vecteur X.

x= a:pas:b ;

Le pas est l’intervalle entre deux points (pas =1 par défaut)

Le vecteur X= [1 2 3 4 5 6 7 8 9 0] peut être défini simplement :

≫ X=1:20 % toutes les valeurs entre 1 et 20 par pas de 1

Une autre alternative, vue précédemment, permet de définir le vecteur en fournissant le nombre de points
voulu.

X=linspace(a,b,n)

Où n est le nombre de points désiré y compris a et b (n=100 par défaut)

Exemple :

≫ X=linspace(1,20,20) % <> X=[1 2 … 20] même vecteur que précédemment

Le vecteur Y est calculé en spécifiant l’expression de la fonction et en prenant en compte des opérateurs
entre vecteur.

Y=f(X) ;

ou f(X) est l’expression de la fonction à tracer.

Ensuite on utilise la fonction plot d’Octave qui réalise le tracé. Dans forme la plus simple, la fonction
plot(x,y) a besoin de deux paramètres d’entrée un vecteur abscisse x, et vecteur ordonné y.

plot (X, Y)

Cette commande joint graphiquement tous les points dont les cordonnées sont contenues dans les vecteurs X
et Y.

Exemple 1 : Tracé de la courbe y = x2 -2x dans l’intervalle [-4, 4]

≫ x = -4:0.1:4; % discrétiser l’axe des x par pas de 0.1

≫ y = x.^2-2*x+1; % calcul des images par la fonction y=x²-2x + 1
≫ plot(x,y) ; % tracé du vecteur y en fonction de x
≫ title('Déplacement vs temps') % Ajouter le titre à la figure
≫ xlabel('Temps’) % Nommer l’axe des abscisses
≫ ylabel('Tension') % Nommer l’axe des ordonnées

Exemple 2 : Tracer la courbe y=sin(x) dans l’intervalle [-π, π]

≫ x=-pi:0.1:pi % discrétiser l’axe des x par pas de 0.1 et

≫ y = sin(x) % calcul du vecteur y
≫ plot(x,y) % tracé de y en fonction de x

10
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Exemple 3 : tracer la fonction 𝑓(𝑥) = 𝑐𝑜𝑠(𝜔1 𝑡) + 𝑐𝑜𝑠(𝜔2 𝑡) , 𝜔1 = 10 , 𝜔2 = 11 𝑟𝑎𝑑. 𝑠 −1 sur

l’intervalle [0, 6𝜋]
≫ w1=10;
≫ w2=11;
≫ t=linspace (0, 6*pi, 500) ;
≫ y = cos(w1*t) + cos(w2*t) ;
≫ plot(t,y) ;

Exercices

Exercice 4 : Tracer l’amplitude 𝑥 𝑒𝑛 𝑓𝑜𝑛𝑐𝑡𝑖𝑜𝑛 𝑑𝑒 𝑡 ∈ [0,20] d’un oscillateur amorti.

𝑥(𝑡) = 𝐴 𝑒 −𝛿𝑡 cos(𝜔𝑎 𝑡 + 𝜙) ; 𝐴 = 10 , 𝜔𝑎 = √𝜔02 − 𝛿 2 𝑎𝑣𝑒𝑐 𝜔0 = 10 , 𝛿 = 1 𝑒𝑡 𝜙 = −0.25

Superposez les graphes, prenez d’autres 𝛿 = 10 𝑒𝑡 𝛿 = 12 par exemple

Exercice 5 : Tracer l’amplitude 𝑒𝑛𝑓𝑜𝑛𝑐𝑡𝑖𝑜𝑛𝑑𝑒 𝜔, 𝜔 ∈ [0, 10]

𝑨𝟎
𝑨(𝝎) = ; 𝑨𝟎 = 𝒆𝟎 𝝎𝟐𝟎 , 𝝎𝟎 = 𝟏𝟎, 𝜹 = 𝟏, 𝒆𝟎 = 𝟎. 𝟒
𝟐
√(𝝎𝟐𝟎 − 𝝎𝟐 ) + 𝟒𝜹𝟐 𝝎𝟐
Conseil :
 Commencer par définir des noms de variables pour les constantes 𝜔, 𝜔0 , 𝑒0 𝑒𝑡 𝛿
 Définir le vecteur des valeurs d’oméga en prenant 100 points dans l’intervalle [0, 10]
 Définir l’expression de A : cette expression est en fonction omega qui est un vecteur. L’expression
sera, donc, appliquée à chaque élément du vecteur omega. Par suite, le résultat sera un vecteur
stocké dans la variable A
 Ayant un vecteur omega qui constituera l’abscisse et le vecteur A constituant l’ordonnée, tracer la
courbe A en fonction omega

Exercice 6 : Reprendre l’exercice 5 avec δ = 0.5 et δ = 2.

Superposer les trois courbes obtenues dans les exercices 5 et 6

Pour superposer les courbes, on utilise la commande :

≫ hold on

Toutes les courbes tracées après cette commande seront superposés sur le même graphe.

Pour arrêter la superposition des courbes :

≫ hold off

11
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Ecriture de scripts simples

Création d’un script :

Lorsqu’on a plusieurs commandes à taper sur la fenêtre de commande pour la réalisation d’une certaine
tâche ; il est judicieux d’utiliser un fichier (script) qui contiendra toutes ces commandes et qui permettrait de
les exécuter en une seule fois. Les commandes seront exécutées dans l’ordre automatiquement.

La création d’un nouveau script se fait sur la fenêtre « Editeur ». Si cette dernière n’est active il faut en co-
chant l’option « Fenêtre-> Afficher l’éditeur » dans le menu principal.

Pour créer un fichier (s’il n’est pas déjà créé par défaut) allez dans le menu de l’éditeur :
Fichier -> Nouveau script
Pour ouvrir un fichier déjà existant Fichier -> ouvrir

Le nom du fichier script doit respecter les règles suivantes :

 Ne doit pas contenir des caractères spéciaux tels que : % / + * , : ; - % …etc.

 Ne doit pas contenir d’espace entre les mots
 Ne doit pas être le nom d’une fonction intégrée dans Octave comme solve, max, min, zeros, …etc.

Exercice 6 : Créer un nouveau script et donner lui un nom « equation2d.m »

Nous allons reprendre l’exemple de la résolution d’une équation de second degré :

𝑥 2 − 13𝑥 + 26 = 0

Dans le script, ajouter les commandes qu’il faut afin de trouver les solutions de l’équation.
(Les commandes que vous avez tapées pour résoudre l’équation dans l’exercice 2)

Pour l’affichage des résultats utilisez la fonction fprintf (format, variables) dont voici la syn-
taxe :
Format : une chaine caractères libre

≫ fprintf(’La solution est x1= %f :’, x1)

Le caractère %f signifie que la variable à afficher à cet endroit est un réel (float en anglais)

%f ou %g : si la variable à afficher est un réel

%d ou %i : si la variable à afficher est un entier
%s : si la variable à afficher est une chaine de caractères

On peut afficher plusieurs variables de types différents :

fprintf(’%f %f \n’, x1, x2 ) ;

Ou de manière plus élaborée

fprintf(’la valeur de x1: %f, la valeur de x2: %f \n’, x1, x2 ) ;

 Il faut qu’il y ait autant de variables que de % dans le format

 Le format est toujours entre deux guillemets ou deux apostrophes

12
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Le caractère \n : pour le retour à la ligne sinon l’affichage continue sur la même ligne à chaque appel à cette
fonction.

Enregistrer le fichier script dans un répertoire sur le disque

Allez dans Fichier->Enregistrer ou cliquez sur l’icône correspondante dans la barre d’outils

Le fichier que vous avez créé aura l’extension ‘.m’ automatiquement

Exécuter le script : allez dans le menu ‘Exécuter le script’ dans le menu de l’éditeur
Ou sur l’icône correspondante dans la barre d’outils de l’éditeur :

A la boite de dialogue suivante, répondre par le bouton du milieu pour rajouter le répertoire du script dans le
chemin de Octave. Le script sera reconnu par Octave et pourra être exécuté depuis la fenêtre de commande
en tapant simplement son nom.

Allez, maintenant, dans fenêtre de commande et l’exécuter sur la ligne de commande en tapant son nom
(sans l’extension)

Modifier les valeurs de : a, b, c dans le script afin de résoudre d’autres cas. Enregistrer puis ré exécuter.

Ecriture de scripts évolués

L’instruction IF

Le mot clé « if » est une instruction de contrôle. Elle permet de faire le test d’une variable logique si elle est
vraie ou fausse. Suivant que le résultat du test est vrai ou faux, elle permet d’exécuter un certain nombre
d’instructions (commandes).

La forme la plus simple du test logique à l’aide de l’instruction « if » est de la forme suivante :

if condition
instruction1
end

exemple : le bloc de code en fond grisé est mettre dans un fichier (nou-
veau script) que vous créez dans Octave.
Ce petit programme demande à l’utilateur d’enter au clavier une valeur de
la variable a. le programme teste ensuite si cette valeur est supérieure
à 0. Si oui, il affiche la phrase ‘a est positif’.

13
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

a=input(‘donner une valeur de a : ’) ; % lire a, à partir du clavier

if (a>0) % test logique
disp(‘a est positif’) ; % instruction (affichage d’une phrase)
end % fin du block if

Si on veut gérer en plus le cas où la condition est fausse

if condition
instruction1
else
instruction2
end

exemple :

a=input(‘donner une valeur de a : ’) ; % lire a à partir du clavier

if (a>0) % test logique
disp(‘a est positif’) ; % instruction1
else
disp(‘a est négatif’) ; % instruction2

end % fin du block if

Il est également possible de rajouter autant de test logique que l’on veut dans un block « if »

if condition1 % si condition1 est vraie, exécuter commande 1

commande 1
elseif condition2 % si condition2 est vraie, exécuter commande 2
commande 2
elseif condition3 % si condition3 est vraie, exécuter commande 3
commande 3

. . . . . . . . . .. . . . . . ..

else % si toutes les conditions sont fausses

commande n % alors exécute la commande n

Exemple :

a=input(‘donner une valeur de a : ’) ; % lire a à partir du clavier

if (a>0) % test logique
disp(‘a est positif’) ; % instruction1
elseif(a<0)
disp(‘a est négatif’) ; % instruction2
else
disp(‘a est nul’) ; % instruction3
end % fin du block if

Exercice 7 : Reprendre le script précédent (resolution d’une equation 2nd degré)

1. Modifier le script pour implémenter un test logique sur le discriminant delta afin de gérer les trois cas
possibles des solutions de l’équation selon le signe de la variable delta
2. Et afficher les deux solutions convenablement à l’aide de la fonction fprintf
3. Exécuter le script en modifiant les constante a, et c afin de vérifier s’il traite tous les cas comme prévu.

14
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Les boucles : For et While

La boucle for

Le mot clé For est utilisée pour répéter un ensemble de commandes un nombre fois connu à l’avance.

La syntaxe de cette structure de contrôle :

for i=first:step:last
command1
command2
command3
command4
. . . . . . . .
. . . . . . . .
end % fin du bloc de commandes à répéter

Les lignes en gras sont forment le bloc de commandes à répéter.

La variable i est: une variable entière : c’est le compteur de répétitions. Elle est incrémentée à chaque exécu-
tion complète du bloc de commandes délimité par les mots clés for et end
first : valeur du début de la variable i
last : la valeur de fin de la variable i
step : le pas d’incrémentation de la variable à chaque exécution du bloc. S’il n’est pas spécifié alors =1

Exemples :
 Afficher les nombres de 1 à 10

≫ for i = 1:10 % le pas n’est pas spécifié donc il est pris égal à 1
≫ disp(i) % affichage du compteur i
≫ end

A chaque répétition affiche le compteur i. Le résultat sera l’affichage tous les nombres entre 1 et 10.

 Afficher les nombres pairs entre 0 et 10 par pas de 2. (les nombres pairs)

for i = 0:2:10
disp(i)
end

Calcul de la somme d’une série de valeurs avec un terme général

Soit par exemple la somme des n premiers entiers

𝑛

𝑆 = ∑𝑖
𝑖=1
La démarche consiste à définir une variable S qui contiendra la somme des n premiers entiers.
Ensuite faire une boucle for entre 1 et n et ajouter la valeur du compteur de répétitions à S et ceci à
chaque répétition.

15
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

S = 0; % initialiser la variable à 0
for i = 1:10 % commence la répétition jusqu’à 10
S = S + i ; % ajouter i à S à chaque répétition
end % fin de la répétition
disp(S) % afficher le résultat

La boucle while

Le mot clé while est utilisé pour répéter l’exécution d’un bloc d’instructions tant qu’une condition
est vraie. Le nombre de répétition n’est pas connu à l’avance.

while expression % expression logique

statements % commande exécuté tant que expression est vraie
end

Attention boucle infinie : en effet, si l’expression logique reste toujours vraie pendant les ré-
pétitions, l’exécution du script ne s’arrêtera pas.
Il n’y a pas de compteur de répétition automatique dans la boucle while. Il doit être implémenté.

Exemples :
 Affichage des nombres inférieurs à 15
≫ i = 0;
≫ while (i< 15 )
≫ i = i + 1;
≫ fprintf('value of a: %d \n', i);
≫ end

 Calcul de la somme des valeurs de 1 à 10.

≫ i = 0; % initiation le compteur à 0
≫ s = 0; % initialiser la somme à 0
≫ while (i <= 10) % tant que le nombre de répétions est inférieur à 15
≫ i = i + 1; % incrémenter le compteur
≫ s = s + i; % à chaque répétition ajouter la valeur de i à s
≫ end

≫ fprintf ('Sum of 1...10 is: %d \n', s);

Exercice 8 : Reprendre le script précédent « equation2d.m »

Modifier le script pour Mettre tout le traitement en boucle

Pour cela, suivre les étapes suivantes :
 Lire les coefficients a, b, c de l’équation 2nd degré depuis le clavier (commande input )
 Faire une boucle while de façon à faire tout le traitement tant que a, b, c sont différent de 0.

Exercices supplémentaires :

16
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Exercice 1 : Calculer la somme de la série suivante avec n=100.

Exercice 2 : Calculer la somme de la série suivante pour n=1000 :

Exercice 3 : Calculer π à l’aide de la formule de Ramanujan (un mathématicien Indien d’exception)

pour n≤10 :

Ecriture de scripts de fonctions

Dans ce chapitre nous allons introduire une technique très importante et très populaire dans la programma-
tion, pour ses grands avantages : les fonctions.

Utilisation d’une fonction :

Une fonction ou une procédure est une portion de code source qui permet de réaliser une certaine tâche (un certain cal-
cul par exemple) bien particulière dont le besoin est récurrent en général.

La fonction reçoit des paramètres d’entrée (dont elle a besoin pour faire sa tache) appelés arguments et retourne le résul-
tat au programme (script) à l’endroit où l’appel à la fonction est fait.

La syntaxe d’appel d’une fonction Matlab :

[r1, r2,...] = functionname(a1,a2,...)

Les arguments : a1, a2 … sont des variables en entrée de la fonction

Les résultats : r1, r2 … sont des variables en sortie dans lesquelles seront stockées les résultats effectués et
retournés donc par la fonction.

Lorsque la fonction retourne une seule valeur on peut écrire (om peut omettre les crochets)
r = functionname(a1,a2,...)

Exemples :
 La fonction sqrt utilise un seule argument X et retourne le résultat dans Y
>> x = 2.3;
>> y=sqrt(x)
y = 1.5166

17
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

 La fonction power (pow2 dans matlab) reçoit deux paramètres en entrée x, y et retourne la valeur 𝑥 𝑦 =
32 qui sera stockée ici, dans la variable z.
>> x=3;
>> y=2;
>> z=power(x,y)
z = 9

Implémentation d’une fonction :

Dans un nouveau fichier, la définition d’une fonction se fait selon la syntaxe suivante :
Le nom du fichier doit être identique au nom de la fonction.

C‘est-à-dire que la fonction funcname sera écrite dans le fichier funcname.m

function [r1, r2,...] = funcname(a1,a2,...)

. . . .
. . . .
end

dans le cas d’une fonction qui retourne une seule valeur :

function r = funcname(a1,a2,...)
. . . .
. . . .
end

Exemple : La fonction qui calcule le carré d’une valeur y = x^2

 On crée un nouveau script dans lequel on écrit la fonction :

function y=carre(x)
y =x^2;
end

Le script est enregistré avec le nom carre dans un le répertoire TP03 par exemple. Le répertoire doit être
dans le chemin de Otave pour se faire il faut aller dans Edit/Path du menu principal puis l’ajouter.

 Dans la fenêtre de commande on calcule le carré de différentes valeurs :

>> carre(3)
ans = 9

 Calcul du carré de la variable a et le stocker dans la variable b

>> a = 3.5000
>> b= carre(a)
b = 12.250

Dans un autre script nommé calcul_carre.m par exemple, on écrit :

x = 2365.2365 ;
y = carre(x) ;
disp(y)

5.5943e+06

Exercice 1 :

18
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Reprendre l’exercice de la résolution de l’équation du second degré.

1. Créer une fonction avec un nom de votre choix.
2. Tester la fonction sur la ligne de commande.
3. Tester la fonction depuis un autre script

Exercice 2 :
On cherche à lisser par une droite les points suivants :
x=[1 2 3 4 5 6 7 8 9 10 11 12] ;
y=[4.3 5.1 5.7 6.3 6.8 7.1 7.2 7.2 7.2 7.2 7.5 7.8 ] ;
Rappel :
Régression linéaire : 𝑦 = 𝑎𝑥 + 𝑏
Avec :
𝑛 𝑛 𝑛
1 𝒄𝒐𝒗𝒙𝒚
𝑥 = ∑ 𝑥𝑖 , 𝑐𝑜𝑣𝑥𝑦 = ∑(𝑥 − 𝑥)(𝑦 − 𝑦) 𝑒𝑡 𝑆𝑥 = ∑(𝑥 − 𝑥)2 → 𝒂 =
2
𝑒𝑡 𝒃 = 𝒚 − 𝒂𝒙
𝑛 𝑺𝟐𝒙
𝑖=1 𝑖=1 𝑖=1
Ecrire une fonction qui calcule les paramètres a et b pour n’importe quel échantillon.
Commencer par écrire un script après l’voir testé convertir en fonction

Conseils :
1- Créer un script nommé regression.m
2- Définir les vecteurs X et Y
3- Définir n=length(X) le nombre d’éléments de X qui est égal au nombre d’éléments de Y aussi.
4- Calculer la moyenne de x qui est la somme divisée par n
5- De même, calculer la moyenne de Y
6- Calculer Sx2 somme des carrés des différences (voir la formule)
7- Calculer la covariance (voir la formule)
8- Calculer a et b. voir formules
9- Tracer Y(X) et y = a x + b, dans le même graphe.

19
USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Partie II :
Travaux pratiques d’analyse
numérique
 Résolutoin d’équations non linéaires
 Interpolation numérique
 Itégration numérique
 Equations différentielle 1er odre
 Equations différentielles 2nd ordre
 Système d’équations linéaires

20
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Résolution des équations non

linéaires du type : 𝒇(𝒙) = 𝟎
Position du problème :
La plupart des équations f(x)=0 rencontrées en physique n’ont pas de solutions analytiques ; d’où le recours à
la résolution numérique.

𝑓(𝑥) : est une fonction à variable réelle. L’idée de la résolution numérique est de construire une suite de va-
leurs numériques (𝑥𝑖 ) qui converge vers une solution r de l’équation. La solution trouvée sera donc une ap-
proximation évaluée avec une certaine tolérance. L’estimation de l’erreur 𝜖 sera prise sous les deux formes
suivantes :

Sur la fonction : |𝑓(𝑥)| ≤ 𝜖 ou sur la solution |𝑥𝑛+1 − 𝑥𝑛 | ≤ 𝜖

Dans ce TP nous allons décrire brièvement les méthodes les plus utilisés et de d’établir un comparatif par
rapport à leur performance.

Préliminaires :
 Afin de comparer leur performance les ; les méthodes seront exécutées avec la même tolérance :
𝜖 = |𝑏 − 𝑎| ≤ 10−6 𝑎𝑢 𝑙𝑖𝑒𝑢 𝑑𝑒 |𝑥𝑛+1 − 𝑥𝑛 | ≤ 10−6 𝑝𝑜𝑢𝑟 𝑠𝑖𝑚𝑝𝑙𝑖𝑐𝑖𝑡é
 Dresser un tableau comparatif comme suit :

La méthode La solution à 10−6 Nombre d’itération Observation

LA dichotomie
Point fixe
Sécante

 On commencera toujours par un tracé de la fonction pour avoir une estimation initiale de sa racine

La méthode de la dichotomie (bissection / bipartition) :

 Condition d’existence de la solution :

Soit 𝑓(𝑥) une fonction définie, continue et strictement monotone sur un intervalle [𝑎, 𝑏].
Si 𝑓(𝑎) × 𝑓(𝑏) < 0 alors 𝑓(𝑥) = 0 possède une solution unique dans [𝑎, 𝑏].
 Convergence : Tout le temps

La méthode numérique :

21
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

1. Avoir une estimation de la solution à l’aide d’un graphique par exemple.

2. Prendre un intervalle quelconque entourant la solution estimée graphiquement (le plus petit entrai-
nera moins d’itérations)
3. Diviser l’intervalle en deux et décider lequel des deux sous intervalles contient la solution. Cette der-
nière sera prise comme étant égale au milieu.
4. Répéter l’opération, l’intervalle se réduit d’une itération à l’autre jusqu’à obtenir la précision voulue.

L’algorithme de la méthode :
Dans un fichier séparé, définir la fonction f1(x)

Initialisation : a, b % les bornes initiales de l’intervalle

Tol % la tolérance = 10^-4
Cpt % compteur des itérations
Traitement : Tant que | b − a| > tol % condition d’arrêt de la boucle
m = (a + b) / 2 % solution = milieu de l’intervalle
Si f1(m)*f1(a) < 0 alors % si m est dans [m, a]
b prend la valeur de m % affecter m à b
Sinon
a prend la valeur de m % affecter m à a
Fin de Si
Incrémenter le compteur Cpt
Fin de Tant que
Sortie : Afficher Cpt, a, b, m

Remarques :
- La valeur absolue est calculée par la fonction matlab abs () : | x| → 𝑎𝑏𝑠 (𝑥)
- Pour simplifier, la tolérance ici porte sur les bornes de l’intervalle et pas sur les mi- intervalles,
ce qui revient au même.

Exercice 1 :
On se propose de calculer le maximum de la densité spectrale du corps noir.

1. Dans un script séparé écrire la fonction f1 : 𝑓1(𝑥) = (3 − 𝑥)𝑒 𝑥 − 3 % max densité d’énergie 𝑢(𝜈, 𝑇)
2. Dans un script séparé écrire la fonction f2 : 𝑓2(𝑥) = (5 − 𝑥)𝑒 𝑥 − 5 % max densité d’énergie 𝑢(𝜆, 𝑇)
3. Dans un autre script programmer la méthode de la dichotomie pour résoudre l’équation 𝑓1(𝑥) = 0
4. Inscrire les informations demandées dans le tableau comparatif

La méthode du point fixe (approximations successives)

Au lieu de chercher la solution de l’équation 𝑓(𝑥) = 0 𝑜𝑛 𝑐ℎ𝑒𝑟𝑐ℎ𝑒 𝑙𝑎 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑑𝑒 𝑔(𝑥) = 𝑥

Il y a de multiple manière de récrire 𝑓(𝑥) = 0 en une telle équation. La plus simple : 𝑔(𝑥) = 𝑥 − 𝑓(𝑥)

 Condition d’existence de la solution :

Si 𝑔 : [𝑎, 𝑏] → [𝑎, 𝑏] dérivable telle que |𝑔′ (𝑥)| < 1 ∀𝑥 ∈ [𝑎, 𝑏]
𝑎𝑙𝑜𝑟𝑠 𝑝𝑜𝑢𝑟 𝑡𝑜𝑢𝑡 𝑥0 ∈ [𝑎, 𝑏], 𝑙𝑎 𝑠𝑢𝑖𝑡𝑒 𝑠𝑢𝑖𝑣𝑎𝑛𝑡𝑒 :
𝑥𝑛+1 = 𝑔(𝑥𝑛 ) 𝑛∈𝑁
𝑐𝑜𝑛𝑣𝑒𝑟𝑔𝑒 𝑣𝑒𝑟𝑠 𝑙 ′ 𝑢𝑛𝑖𝑞𝑢𝑒 𝑝𝑜𝑖𝑛𝑡 𝑓𝑖𝑥𝑒 𝑑𝑒 𝑔.

22
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Convergence : |𝑔′ (𝑥)| < 1 ∀𝑥 ∈ [𝑎, 𝑏]

Il est à noter que pour la même équation initiale, certaines réécritures en une équation 𝑔(𝑥) = 𝑥
peuvent converger et que d’autres ne convergent pas

L’algorithme de la méthode :

Initialisation : % une première estimation de la solution

𝑥0
tol
% la tolérance
Cpt % compteur des itérations
𝑥1 < −𝑔(𝑥0 ) % 1er calcul de x1 pour calculer l’erreur
Traitement : Tant que | x1 − x0 | > tol % condition d’arrêt de la boucle
𝑥0 < −𝑥1 % x0 devient x1 pour le calcul suivant
𝑥1 < −𝑔(𝑥0 ) % refaire l’opération
Cpt <- Cpt+1 % Incrémenter le compteur Cpt
Fin de Tant que
Sortie : Afficher Cpt, x1

Conseil : Tracer la fonction g(x) en superposition avec la droite y=x pour une estimation de la solu-
tion.

Exercice 2 :
1. Créer une fonction g définie comme suit : 𝑔(𝑥) = 3(1 − 𝑒 −𝑥 )
2. Tracer sur la même figure la fonction g et la droite 𝑦 = 𝑥
3. Vérifier la convergence de g et prendre une estimation du départ.
4. Ecrire un nouveau script nommé point_fixe.m et programmer la méthode
5. Exécuter le script pour la fonction f1 et inscrire les informations demandées.

La méthode de Newton (Newton-Raphson)

 Condition d’existence de la solution :

𝑓 ∶ 𝑓𝑜𝑛𝑐𝑡𝑖𝑜𝑛 𝑑é𝑓𝑖𝑛𝑖𝑒 𝑒𝑡 𝑑𝑒𝑢𝑥 𝑓𝑜𝑖𝑠 𝑑é𝑟𝑖𝑣𝑎𝑏𝑙𝑒 𝑠𝑢𝑟 𝑢𝑛 𝑖𝑛𝑡𝑒𝑟𝑣𝑎𝑙𝑙𝑒 𝐼 = [𝑎, 𝑏], et qu’elle possède une racine
r dans I

𝑓(𝑥𝑛 )
𝑥𝑛+1 = 𝑥𝑛 − 𝑛∈𝑁
𝑓 ′ (𝑥𝑛 )

Convergence : f définie et deux fois dérivable tel que 𝑓 ′ (𝑥) ≠ 0 au voisinage de r

L’algorithme de la méthode :

Définir la fonction qui calcul les fonctions 𝑓(𝑥) 𝑒𝑡 𝑓 ′ (𝑥)

Initialisation : 𝑥0 , 𝑡𝑜𝑙, 𝑐𝑝𝑡
Calculer x1 = x0 − f(x0 )⁄f ′ (𝑥0 )
Traitement : Tant que |x1 − x0 | > tol
Affecter x1 à x0

23
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Calculer x1 = x0 − f(x0 )⁄f ′ (x0 )

Affecter x1 à x0
Incrémenter cpt
Fin Tant que
Afficher cpt, xn

La méthode de la sécante

C’est la méthode de Newton modifié. La dérivée est intégrée dans le schéma numérique.
𝑥𝑛 − 𝑥𝑛−1
𝑥𝑛+1 = 𝑥𝑛 − 𝑓(𝑥𝑛 )
𝑓(𝑥𝑛 ) − 𝑓(𝑥𝑛−1 )

Algorithme de la méthode :
Initialisation : 𝑥0 = 𝑎, 𝑥1 = 𝑏, 𝑡𝑜𝑙 = 10−6 , 𝑐𝑝𝑡 = 0
Traitement : 𝑇𝐴𝑁𝑇 𝑄𝑈𝐸 |𝑥1 − 𝑥0 | > 𝑒𝑟𝑟
𝑥1 −𝑥0
𝐶𝑎𝑙𝑐𝑢𝑙𝑒𝑟 𝑥2 = 𝑥1 − 𝑓(𝑥1 ) 𝑓(𝑥
1 )−𝑓(𝑥0 )

𝐴𝑓𝑓𝑒𝑐𝑡𝑒𝑟 𝑥1 à 𝑥0
𝐴𝑓𝑓𝑒𝑐𝑡𝑒𝑟 𝑥2 à 𝑥1
𝐼𝑛𝑐𝑟é𝑚𝑒𝑛𝑡𝑒𝑟 𝑙𝑒 𝑐𝑜𝑚𝑝𝑡𝑒𝑢𝑟 𝑐𝑝𝑡
𝐹𝐼𝑁 𝑇𝐴𝑁𝑇 𝑄𝑈𝐸
𝐴𝑓𝑓𝑖𝑐ℎ𝑒𝑟 𝑐𝑝𝑡, 𝑥𝑛

Exercice 3:
1. Créer un nouveau script nommé secante.m, dans lequel vous programmez la méthode de la sécante
2. Exécuter la méthode pour la fonction f1.
3. Finir le comparatif

Exercice 4:
1. Convertir le script en une fonction Octave qui reçoit en entrée : le nom de la fonction, la tolérance
et éventuellement une première estimation 𝑥0 .
2. Afficher le résultat : le nombre d’itération et la solution de l’équation.
3. Exécuter la fonction pour les fonctions suivantes :
𝑓2 (𝑥) = 𝐽1 (𝑥) dans l intervalle [0, 5] ou 𝐽: la fonction de Bessel (𝑏𝑒𝑠𝑠𝑒𝑙𝑗(0, 𝑥) )
4. Tracer les fonctions de Bessel : 𝐽1 , 𝐽2 , 𝐽3 , 𝐽4 , 𝐽5 sur l’intervalle [0, 20] pour l’estimation des ra-
cines
5. Prendre une 1ere estimation des racines et trouver ces premières racines avec une tolérance 𝑡𝑜𝑙 =
10−4 de ces fonctions en utilisant votre fonction Octave : secante.

24
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Interpolation polynomiale
Introduction

L’interpolation est une approximation d’une fonction inconnue ou non à partir d’une série de points
qui sont des mesures de cette fonction. Autrement dit, représenter cette fonction au mieux à l’aide
d’un polynôme.
Le but de l’interpolation est d’avoir un moyen pour prédire les valeurs de cette fonction en des
points ou on ne dispose pas de points de mesures. Les points interpolés doivent être à l’intérieur de
l’intervalle des points mesurés.
Deux méthodes sont présentées : la méthode de Lagrange et la méthode de Newton.

Interpolation de Lagrange
Introduction
Approximation d’une fonction f sur un intervalle [a, b] ou on dispose de n points (xi, yi) par un polynôme de
degré n tel que 𝑝(𝑥𝑖 ) = 𝑓(𝑥𝑖 ) , 𝑖 = 1, … 𝑛 tel que 𝑎 ≤ 𝑥1 ≤ 𝑥2 ≤ ⋯ ≤ 𝑥𝑛

Il existe un polynôme unique qui vérifie ce critère à savoir : le polynôme d’interpolation de Lagrange :
𝑛

𝑃(𝑥) = ∑ 𝑓(𝑥𝑘 )𝐿𝑘 (𝑥)

𝑘=1

Ou 𝐿𝑘 (𝑥) 𝑘 = 1, … 𝑛 est la base des polynômes de Lagrange :

𝑛
𝑥 − 𝑥𝑖
𝐿𝑘 (𝑥) = ∏
𝑥𝑘 − 𝑥𝑖
𝑖=1,𝑖≠𝑘

On constate que la modification, l’ajout ou la suppression d’un point quelconque nécessite de recal-
culer de nouveau tous les polynômes de la base.

Exemple :
Soit les mesures suivantes :

X 1 2 3
Y 0.3 0.5 0.2

Cherchons l’interpolation de Lagrange en 𝑥𝑝 = 2.5

𝑥𝑝 − 𝑥2 𝑥𝑝 − 𝑥3 2.5 − 2 2.5 − 3
𝐿1 (𝑥𝑝 ) = × = × = −0.125
𝑥1 − 𝑥2 𝑥1 − 𝑥3 1−2 1−3

25
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

𝑥𝑝 − 𝑥1 𝑥𝑝 − 𝑥3 2.5 − 1 2.5 − 3
𝐿2 (𝑥𝑝 ) = × = × = 0.75
𝑥2 − 𝑥1 𝑥2 − 𝑥3 2−1 2−3

𝑥𝑝 − 𝑥1 𝑥𝑝 − 𝑥2 2.5 − 1 2.5 − 2
𝐿3 (𝑥𝑝 ) = × = × = 0.375
𝑥3 − 𝑥1 𝑥3 − 𝑥2 3−1 3−2

𝑃(𝑥) = 𝑓(𝑥1 )𝐿1 (𝑥𝑝 ) + 𝑓(𝑥2 )𝐿2 (𝑥𝑝 ) + 𝑓(𝑥3 )𝐿3 (𝑥𝑝 ) (1)

= 0.3 × −0.125 + 0.5 × 0.75 + 0.2 × 0.375 = 0.4125

Pour une valeur de x donnée 𝑥 = 𝑥𝑝 , la méthode numérique sous Octave / Matlab consiste à

1- Calculer les valeurs 𝐿1 (𝑥𝑝 ) , 𝐿2 (𝑥𝑝 ) , 𝐿3 (𝑥𝑝 ) , … 𝐿𝑛 (𝑥𝑝 ) qui seront stockés dans un vecteur qu’on
appellera L : 𝐿 = [𝐿1 (𝑥𝑝 ), 𝐿2 (𝑥𝑝 ), … , 𝐿𝑛 (𝑥𝑝 )]

2- De façon similaire à (1) on obtient 𝑝(𝑥𝑝 ) = 𝑓(𝑥1 )𝐿1 (𝑥) + 𝑓(𝑥2 )𝐿2 (𝑥) + ⋯ + 𝑓(𝑥𝑛 )𝐿𝑛 (𝑥) qui n’est
rien d’autre que le produit des deux vecteurs L par le vecteur f avec f = [𝑓(𝑥1 ), 𝑓(𝑥2 ), … , 𝑓(𝑥𝑛 )]

Algorithme
Données :
X : vecteur contenant les abscisses des points de mesure (𝑙𝑒𝑠 𝑥𝑖 , 𝑖 = 1 … 𝑛)
Y : vecteur contenant les ordonnées des points de mesure (𝑙𝑒𝑠 𝑓(𝑥𝑖 ) , 𝑖 = 1 … 𝑛)
xp : La valeur test de l’interpolation

% calcul des polynômes L(k) k=1...n

Pour k de 1 à n
L(k) <- 1
Pour i de 1 à n
Si k != i
L(k) <- L(k) * (x - X(i)) / (X(k) - X(i));
Fin Si
Fin Pour
Fin Pour
% calcul de l’interpolation de Lagrange en x=xp
P <- somme (Y.*L) ; % la variable P contient le résultat

Remarque : La fonction sum(X) retourne la somme des éléments du vecteur X.

Travail demandé
Les mesures de l’amplitude de la tension aux bornes du condensateur d’un oscillateur RLC en régime forcé
sont inscrites dans le tableau suivant :

f(kHz) 5.00 6.00 7.00 7.50 8.00 8.50 9.00 9.5 10 10.5 11 12 13
Vc(V) 0.93 1.24 1.99 3.01 6.21 9.43 3.52 2.01 1.37 1.03 0.81 0.56 0.42

On souhaite avoir une meilleure estimation de la fréquence de résonnance.

26
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

1- Tracer les points Vc en fonction de f et observez l’allure de la courbe notamment à l’alentour de la

résonnance.
2- Ecrire un script Matlab du nom interp_lagrange.m pour évaluer la tension aux fréquences f=8.1
3- Convertir le script en une fonction qui accepte deux vecteurs en entrée : X, Y et une valeur d’interpo-
lation xp. La fonction doit donner en retour la valeur d’interpolation (tension) en ce point.
4- Calculer l’interpolation pour f= 8.2, 8.3 et 8.4. Conclure ?

Interpolation de Newton
Introduction
De la même façon que l’interpolation de Lagrange, l’interpolation de Newton approche une fonction f sur un
intervalle [a, b] ou on dispose de n points (xi, yi) par un polynôme de degré n tel que 𝑝(𝑥𝑖 ) = 𝑓(𝑥𝑖 ) , 𝑖 =
1, … 𝑛 tel que 𝑎 ≤ 𝑥1 ≤ 𝑥2 ≤ ⋯ ≤ 𝑥𝑛

De même un polynôme unique existe et qui sera identique à celui de Lagrange. La différence est dans la ma-
nière d’obtenir ce polynôme.
𝑛

𝑃𝑛 (𝑥) = ∑ 𝛼𝑘 𝑒𝑘 (𝑥)
𝑘=1

 𝑒𝑘 (𝑥) : est le polynôme de la base de Newton d’ordre k.

𝑘−1

𝑒𝑘 (𝑥) = ∏(𝑥 − 𝑥𝑖 ) , 𝑘 = 1, … 𝑛
𝑖=1

 Les coefficients sont obtenus par un calcul des différences divisées : 𝛼𝑘 , 𝑘 = 1, … 𝑛

Le polynôme d’interpolation de Newton s’écrit alors
𝑛

𝑃𝑛 (𝑥) = ∑ 𝛼𝑘 𝑒𝑘 (𝑥)
𝑘=1

Les différences divisées

Newton définit les différences divisées comme suit :
Etant données n points de mesures : (xi , 𝑓𝑖 ) 𝑖 = 1, … 𝑛

𝑓2 − 𝑓1 𝑓3 − 𝑓2 𝑓𝑛 − 𝑓𝑛−1
𝑓21 = , 𝑓32 = , … 𝑓𝑛(𝑛−1) =
𝑥2 − 𝑥1 𝑥3 − 𝑥2 𝑥𝑛 − 𝑥𝑛−1

𝑓32 − 𝑓21 𝑓43 − 𝑓32 𝑓𝑛(𝑛−1) − 𝑓(𝑛−1)(𝑛−2)

𝑓321 = , 𝑓432 = , … 𝑓𝑛(𝑛−1)(𝑛−2) =
𝑥3 − 𝑥1 𝑥4 − 𝑥2 𝑥𝑛 − 𝑥𝑛−2

𝑓432 − 𝑓321 𝑓543 − 𝑓432

𝑓4321 = , 𝑓5432 = , … 𝑒𝑡𝑐
𝑥4 − 𝑥1 𝑥5 − 𝑥2
…………

27
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Les premiers éléments de chaque série des différences divisées constituent les coefficients αk du polynôme
d’interpolation
𝑓𝑘….2 − 𝑓𝑘−1….1
𝑓𝑘….1 = 𝑎𝑣𝑒𝑐 𝑖 = 1 … 𝑛
𝑥𝑘 − 𝑥1

Le calcul de ces éléments (diagonaux) nécessite le calcul de tous les autres différences divisées

Exemple : calcul des différences divisées inférieures pour n=7

𝑥1 𝑓1

𝑥2 𝑓2 𝑓2 − 𝑓1
𝑓21 =
𝑥2 − 𝑥1
𝑥3 𝑓3 𝑓3 − 𝑓2 𝑓32 − 𝑓21
𝑓32 = 𝑓321 =
𝑥3 − 𝑥2 𝑥3 − 𝑥1
𝑥4 𝑓4 𝑓43 𝑓43 − 𝑓32 𝑓4321
𝑓432 =
𝑥4 − 𝑥2
𝑥5 𝑓5 𝑓54 𝑓543 𝑓5432 𝑓54321

𝑥6 𝑓6 𝑓65 𝑓654 𝑓6543 𝑓65432 𝑓654321

𝑥7 𝑓7 𝑓76 𝑓765 𝑓7654 𝑓76543 𝑓765432 𝑓7654321

Les coefficients 𝛼𝑘 sont : 𝛼1 = 𝑓1 , 𝛼2 = 𝑓21 , 𝛼3 = 𝑓321 … 𝛼7 = 𝑓7654321

Traitement numérique
Pour calculer les coefficients du polynôme d’interpolation, il est pratique, de constituer une table, comme
dans l’exemple précédent, de toutes les différences divisées. Ensuite, prendre les éléments de la diagonale
comme étant les coefficients recherchés.

Pour le calcul des différences divisées sous Matlab/Octave, nous avons besoin d’un terme général en fonc-
tion de i,j les numéros de ligne et de colonne de la matrice des différences divisées:

i 1 2 3 4 5 6 7

j
1 𝑑1 = 𝑓1

2 𝑑2 = 𝑓2 𝑓2 − 𝑓1
𝑑22 =
𝑥2 − 𝑥1
3 𝑑3 = 𝑓3 𝑓3 − 𝑓2 𝑑32 − 𝑑22
𝑑32 = 𝑑33 =
𝑥3 − 𝑥2 𝑥3 − 𝑥1
4 𝑑4 = 𝑓4 𝑑42 𝑑42 − 𝑑32 𝑑44
𝑑43 =
𝑥4 − 𝑥2
5 𝑑5 = 𝑓5 𝑑52 𝑑53 𝑑54 𝑑54 − 𝑑44
𝑑55 =
𝑥5 − 𝑥1

28
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

6 𝑑6 = 𝑓6 𝑑62 𝑑63 𝑑64 𝑑64 − 𝑑54 𝑑66

𝑑65 =
𝑥6 − 𝑥2
7 𝑑7 = 𝑓7 𝑑72 𝑑73 𝑑74 𝑑74 − 𝑑64 𝑑76 𝑑77
𝑑75 =
𝑥7 − 𝑥3

Illustration du calcul des éléments de la matrice des différences

Le terme général de l’élément 𝑑𝑖𝑗 peut s’écrire sous la forme suivante :

𝑑𝑖,𝑗−1 − 𝑑𝑖−1,𝑗−1
𝑑𝑖𝑗 = 𝑖 = 2 … 𝑛, 𝑗 = 2 … 𝑛
𝑥𝑖 − 𝑥𝑖−𝑗+1
Les éléments 𝑑𝑖𝑗 sont les éléments de la matrice des différences divisées D.

Le schéma de Horner
Le polynôme
𝑛

𝑃𝑛 (𝑥) = ∑ 𝛼𝑘 𝑒𝑘 (𝑥)
𝑘=1
Peut s’écrire sous la forme :

𝑃𝑛 (𝑥) = 𝛼1 + 𝛼2 (𝑥 − 𝑥1 ) + 𝛼3 (𝑥 − 𝑥1 )(𝑥 − 𝑥2 ) + 𝛼4 (𝑥 − 𝑥1 )(𝑥 − 𝑥2 )(𝑥 − 𝑥3 ) … 𝛼𝑛 (𝑥 − 𝑥1 ) … (𝑥 − 𝑥𝑛−1 )

𝑃𝑛 (𝑥) = 𝛼1 + (𝑥 − 𝑥1 ) (𝛼2 + (𝑥 − 𝑥2 ) (𝛼3 + (𝑥 − 𝑥3 )(𝛼4 … 𝑎𝑛−2 + (𝑥 − 𝑥𝑛−2 )(𝛼𝑛−1 + (𝑥 − 𝑥𝑛−1 )𝛼𝑛 ))))

En faisant le produit à partir de de l’indice le plus grand :

l’opération × (𝑥 − 𝑥𝑘 ) + 𝛼𝑘 est récurrente sauf pour le dernier terme 𝛼𝑛 .
Par suite, si la valeur initiale du polynôme est égale à 𝛼𝑛 et qu’en faisant varier k de n-1 à 1 on évalue le polynôme

Intérêt : Une seule boucle permet d’évaluer 𝑃𝑛 (𝑥) au lieu de deux et donc nombre d’opéra-
tion optimale.

Exemple :
𝑃𝑛 (𝑥) = 𝑎𝑥 3 + 𝑏𝑥 2 + 𝑐𝑥 + 𝑑 → 𝑛𝑜𝑚𝑏𝑟𝑒 𝑑′ 𝑜𝑝é𝑟𝑎𝑡𝑖𝑜𝑛 = 9

𝑃𝑛 (𝑥) = 𝑑 + 𝑥(𝑐 + 𝑥(𝑏 + 𝑎𝑥)) → 𝑛𝑜𝑚𝑏𝑟𝑒 𝑑′ 𝑜𝑝é𝑟𝑎𝑡𝑖𝑜𝑛𝑠 = 6

29
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Pour évaluer le polynôme de Newton nous allons opter pour le schéma de Horner pour son intérêt en termes
de performance.

Algorithme
Données : vecteur X, Y contenant les points de mesures
Xp : abscisse pour lequel l’interpolation est évaluée

1. Définir une matrice D remplie de zéros initialement.

2. Affecter les Yi la première colonne : D (:1) =Y
3. Faire une première boucle sur les lignes i allant de 2 à n
4. Faire une deuxième boucle sur les colonnes j allant de 2 à n
5. Calculer les éléments de la matrice D(i,j)= . . .
6. Après les deux boucles, le vecteur a, sera égal à la diagonale de
D.
7. Calculer l’interpolation à l’aide du schéma de Horner suivant :

% P est la variable qui contiendra l’interpolation

P = a(N); % initialement P=a(n)
for k=N-1:-1:1 % de n-1 à 1 par pas de -1
P = P*(x - X(k)) + a(k);
end
Afficher P
Travail demandé
Répondre aux mêmes questions que celles de la partie « interpolation de Lagrange »

30
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Intégration numérique
Introduction

Deux cas, au mois, nécessitent le recours aux méthodes numériques à la recherche de la meilleure approxi-
mation de l’intégrale :
𝑏
∫ 𝑓(𝑥)
𝑎
𝜋
 La fonction 𝑓(𝑥) n’a pas de primitive analytique. Exemple : ∫0 cos 𝑥 2 𝑑𝑥
 La fonction est inconnue mais on en dispose d’une série de mesure.
 Recherche d’une moyenne d’un signal représentant une série de mesure :

𝑏
1
𝑓̅ = ∫ 𝑓(𝑥)
𝑏−𝑎 𝑎

Deux méthodes seront présentées dans ce chapitre à savoir : la méthode des trapèzes et la méthode de Simp-
son.

La méthode des trapèzes

La méthode du trapèze consiste à approximer

l’air sous la courbe par la surface du trapèze dont
la petite base est égale à f(a) et la grande base est
égale f(b).

𝑏
(𝑓(𝑎) + 𝑓(𝑏)) × (𝑏 − 𝑎)
𝐼 = ∫ 𝑓(𝑥)𝑑𝑥 ≅
𝑎 2

La méthode numérique consiste, donc, à subdiviser l’intervalle [a, b] en de petits sous-intervalles et calculer,
ensuite, l’air de chaque sous intervalle par la formule du trapèze ci-dessus.

31
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Le cumul de toutes ces airs, donne une approxi-

mation de l’intégrale totale

𝑛
ℎ
𝐼 = ∑[𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑖+1 )]
2
𝑖=1

𝑏−𝑎
ℎ=
𝑛

L’algorithme
L’intégrale peut s’écrire sous la forme suivante :
𝑛 𝑛 𝑛+1
ℎ ℎ
𝐼 = ∑[𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑖+1 )] = 𝐼 = (∑𝑓(𝑥𝑖 ) + ∑ 𝑓(𝑥𝑖 ))
2 2
𝑖=1 𝑖=1 𝑖=2

Données :
a, b % Les bornes de l’intégrale
n % nombre de sous intervalles
f(x) % la définition de la fonction à intégrer

1. Définir les bornes de l’intervalle d’intégration : a et b

2. Définir la fonction à intégrer ex : f=@(x) expression ;
3. Calculer n le nombre des sous intervalles
4. Définir un vecteur X des abscisses des n+1 points
5. Définir un vecteur Y des images des éléments de X par la fonction f
6. A l’aide de la fonction SUM de Matlab / Octave calculer les deux
sommations et en déduire la valeur de l’intégrale.

Attention :
𝑛

∑𝑓(𝑥𝑖 ) Correspond à sum(Y)-f(b)

𝑖=1
𝑛+1

∑𝑓(𝑥𝑖 ) Correspond à sum(Y)-f(a)

𝑖=2

Travail demandé :

Créer un script dans lequel il faut implémenter la méthode des trapèzes.

Le tester pour l’intégrale I1.
Convertir le script en une fonction matlab appelée trapaze qui reçoit les paramètres f, a, b, n par argument.

Tester la fonction comme suit :

Sur la fenêtre de commande definissez a, b, n et
f = @(x) . . .

32
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Ensuite calculer les intégrales suivantes à l’aide de cette fonction matlab.

2 𝜋
𝐼1 = ∫ √𝑥𝑑𝑥 𝑛 = 10 𝐼3 = ∫ cos2 𝑥 𝑑𝑥 𝑛 = 10
1 0
2 𝜋
𝐼2 = ∫ √1 + 𝑒 𝑥 𝑑𝑥 𝑛 = 20 𝐼4 = ∫ 𝑐𝑜𝑠 𝑥 2 𝑑𝑥 𝑛 = 10
0 0

Test d’autres valeurs pour n et remarquer la constater la précision du calcul

33
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

La méthode de Simpson
La méthode de Simpson consiste à approximer

La fonction f par un polynôme d’interpolation

par la méthode de Lagrange ou de Newton puis
effectuer l’intégration analytique.

Interpolation de Lagrange d’ordre 2 sur [a, b]:

(𝑥 − 𝑚)(𝑥 − 𝑏)
𝑃(𝑥) = 𝑓(𝑚)
(𝑎 − 𝑚)(𝑎 − 𝑏)
(𝑥 − 𝑏)(𝑥 − 𝑎)
+ 𝑓(𝑚)
(𝑚 − 𝑏)(𝑚 − 𝑎)
(𝑥 − 𝑚)(𝑥 − 𝑎)
+ 𝑓(𝑏)
(𝑏 − 𝑚)(𝑏 − 𝑎)

L’intégrale vaut :
𝑏 𝑏
𝑏−𝑎
∫ 𝑓(𝑥)𝑑𝑥 ≅ ∫ 𝑃(𝑥)𝑑𝑥 = [𝑓(𝑎) + 4𝑓(𝑚) + 𝑓(𝑏)]
𝑎 𝑎 6

La méthode numérique consiste donc à subdiviser l’intervalle [a, b] en n sous intervalles :

𝑥𝑖 = 𝑎 + (𝑖 − 1)ℎ 𝑖 = 1 … 𝑛

De calculer, ensuite, l’intégrale de la fonction dans chaque sous intervalles à l’aide de la formule de Simpson
ci-dessus, ceci donnera :

𝑥𝑖+2 𝑥𝑖+2
ℎ 𝑥𝑖 + 𝑥𝑖+1
∫ 𝑓(𝑥)𝑑𝑥 ≅ ∫ 𝑃(𝑥)𝑑𝑥 = [𝑓(𝑥𝑖 ) + 4𝑓(𝑚𝑖 ) + 𝑓(𝑥𝑖+1 )] , 𝑚𝑖 =
𝑥𝑖 𝑥𝑖 6 2

Le cumul sur tous les intervalles donnera l’évaluation de l’intégrale de la fonction f sur l’intervalle [a, b]

𝑏 𝑛
ℎ
∫ 𝑓(𝑥)𝑑𝑥 ≅ ∑ [𝑓(𝑥𝑖 ) + 4𝑓(𝑚𝑖 ) + 𝑓(𝑥𝑖+1 )]
𝑎 6
𝑖=1

L’algorithme
L’intégrale peut s’écrire sous la forme :
𝑛−1 𝑛 𝑛 𝑛+1
ℎ ℎ
𝐼 = ∑ [𝑓(𝑥𝑖 ) + 4𝑓(𝑚𝑖 ) + 𝑓(𝑥𝑖+1 )] = (∑ 𝑓(𝑥𝑖 ) + 4 ∑ 𝑓(𝑚𝑖 ) + ∑ 𝑓(𝑥𝑖 ))
6 6
𝑖=1 𝑖=1 𝑖=1 𝑖=2
𝑝𝑎𝑠=2

Pour implémenter la méthode de Simpson suivre les étapes suivantes :

34
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Données :
a, b % Les bornes de l’intégrale
n % nombre de sous intervalles donc n+1 points
f(x) % la définition de la fonction à intégrer

1. Définir les bornes de l’intervalle d’intégration : a et b

2. Définir la fonction à intégrer ex : f=@(x) expression ;
3. Calculer n le nombre des sous intervalles
1. Définir un vecteur X des abscisses des n+1 points
2. Définir un vecteur Y des images des éléments de X
3. Définir le vecteur M des points milieux
4. Définir le vecteur YM des images des éléments de M
5. A l’aide de la fonction SUM de Matlab / Octave calculer les diffé-
rentes sommations et en déduire la valeur de l’intégrale.

Travail demandé :
Créer une fonction matlab appelée simpson qui reçoit les paramètres f, a, b, n par argument.

A l’aide de la fonction : I=simpson(f,a,b,n)

Calculer les mêmes intégrales de l’exercice 1 et observer.

Exercice d’application

Modifier la(les) fonction(s) pour prendre en argument un ensemble de points :

I=simpson (X, Y) ou I=trapeze(X, Y)

Puis, répondre à l’exercice suivant :

On lance une fusée verticalement du sol et l’on mesure pendant 80 secondes l’accélération γ, les mesures
sont inscrites dans le tableau suivant :

t(s) 0 10 20 30 40 50 60 70 80
γ (m.s-2 30 31.63 33.44 35.47 37.75 40.33 43.29 46.70 50.67

Exprimer la vitesse de la fusée à l’instant t=80s, puis la calculer à l’aide des fonctions que vous avez pro-
grammé.

35
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Equations différentielles de premier ordre

La résolution d’équations différentielles analytiquement est une exception. En effet, pour la quasi majorité des
Equations différentielles, on peut démontrer l’existence de la solution mais sans pouvoir déterminer cette so-
lution analytiquement, d’où l’appel à la résolution numérique.

Deux sortes de problèmes se posent :

 On dispose d’une condition initiale (𝑡0 , 𝑦0 ) et d’une expression 𝑦 ′ = 𝑓(𝑡, 𝑦)
 On dispose d’une série de mesure de (𝑡𝑖 , 𝑦𝑖 ) 𝑖 = 1 … . 𝑛

Nous allons présenter dans ce qui suit les méthodes suivantes :

 La méthode d’Euler
 La méthode de Runge-Kutta d’ordre 1
 La méthode de Runge-Kutta d’ordre 2

La méthode générale consiste à écrire l’équation différentielle sous la forme de Cauchy :

𝑦 ′ = 𝑓(𝑡, 𝑦(𝑡))
{ ∀ 𝑡 𝜖 [𝑡0 , 𝑡0 + 𝑇] (1)
𝑦(𝑡0 ) = 𝑦0

La méthode d’Euler
Cette méthode est loin d’être la plus performante, elle est présentée pour une raison didactique par simplicité.

Soit la forme de Cauchy (1), la méthode consiste à user du développement limité de la fonction y :

𝒅𝒚
𝑦(𝑡 + Δ𝑡) = 𝑦(𝑡) + Δ𝑡 + 𝑂(Δ𝑡 2 )
𝒅𝒕
𝒅𝒚
Utilisons l’expression de 𝒅𝒕 définie dans (1) on obtient :

𝑦(𝑡 + Δ𝑡) = 𝑦(𝑡) + Δ𝑡𝒇(𝒕, 𝒚) (2)

Le schéma peut être vu comme une intégration numérique par la méthode des
rectangles. En effet le terme Δ𝑡𝒇(𝒕, 𝒚) représente l’air sous la courbe f dans
l’intervalle Δ𝑡.

L’intégration de f qui est la solution du problème est le cumul des intégrales de f dans les sous-intervalles Δ𝑡.

Le schéma numérique
Si on discrétise l’axe de la variable t :

𝑡 = 𝑡0 + 𝑖Δ𝑡 𝑖 = 1, 𝑛 𝑛 = 𝑇/Δ𝑡

Ou T est l’intervalle d’intégration de l’équation différentielle.

Si on applique le schéma (2) pour calculer (i+1)eme intervalle connaissant le ieme, on obtient :

36
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

𝑦𝑖+1 = 𝑦𝑖 + Δ𝑡 𝒇(𝒕𝒊 , 𝒚𝒊 ) , 𝑖 = 1 … 𝑛 (3)

La valeur y1 est donné dans la formulation de Cauchy.

L’algorithme
Pour implémenter la méthode d’Euler suivre les étapes suivantes :
Il est usage dans la résolution numérique des équations différentielles de commencer par définir un pas de
temps puis en déduire le nombre de points à utiliser au lieu de l’inverse (comme dans le cas de l’interpolation
ou l’intégration numérique). En effet, dans la résolution des problèmes de la physique, un présentiment sur
l’ordre de grandeur du pas du temps est primordial.

Commencer par définir les données du problème :

1. Définir un intervalle d’intégration c’est-à-dire : 𝑡0 𝑒𝑡 𝑡𝑚𝑎𝑥
2. Définir le pas d’intégration dt

Selon la forme de Cauchy il faut définir une fonction à deux variable yp(t,y) pour chaque problème à traiter.
Cette fonction peut accepter indifféremment des scalaires ou des vecteurs.
3. Définir la fonction yp

Comme pour la plupart des méthodes numériques, il faut commencer par discrétiser la variable d’intégration.
4. Définir un vecteur t qui va de t0 à tmax par pas de dt.

Nous avons besoin du nombre de points pris en considération pour réaliser la boucle d’intégration.
5. Définir n le nombre d’éléments du vecteur t.

Les valeurs recherchées (solution de l’équation différentielle à résoudre) sont les y(i). Le premier élé-
ment du vecteur y est y0.
6. Définir le premier élément du vecteur y.

Maintenant, il faut réaliser la boucle pour calculer tous les éléments du vecteur y dans l’ordre : y(2),
y(3), …y(n) selon le schéma numérique (3).

Trouver les bornes de la boucle pour ce schéma puis l’implémenter

7. Réaliser la boucle d’intégration

Après la terminaison de la boucle, nous avons besoin de voir l’allure de la solution y.

8. Tracer y(t)

Travail demandé
Créer un script dans Octave qui implémente la méthode d’Euler en suivant les étapes de l’algorithme ci-des-
sus pour résoudre numériquement l’équation différentielle suivante :

𝑦 ′ = −𝑡𝑦 𝑡 ∈ [0, 10]

{
𝑦(0) = 1

1- Prenez un pas d’intégration égal à 0.1

2- Tracer la courbe de la solution y(t)
3- Superposez la courbe de la solution exacte que vous déterminer analytiquement
4- Superposez les deux courbes pour comparer
5- Spécifiez les données du problème, ensuite convertir le script en un script d’une fonction euler(…) .
Cette fonction reçoit ces données en paramètres d’entrée et retourne un vecteur y contenant la solu-
tion de l’équation différentielle.

37
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Exécuter la fonction pour le problème suivant :

𝑦 ′ = 𝑦 + 𝑒 2𝑡 𝑡 ∈ [0, 10]
{
𝑦(0) = 1

La méthode de Runge-Kutta ordre 2

C’est une amélioration de la méthode d’Euler. Elle consiste à

estimer l’intégrale de f par un trapèze au lieu du rectangle.
𝑏 𝑏
′
𝑏−𝑎
∫ 𝑦 𝑑𝑡 = ∫ 𝑓(𝑡, 𝑦)𝑑𝑡 = [𝑓(𝑎) + 𝑓(𝑏)]
𝑎 𝑎 2

Pour l’intervalle i :
𝑡𝑛+1
ℎ
𝑦𝑛+1 − 𝑦𝑛 = ∫ 𝑓(𝑡, 𝑦)𝑑𝑡 = [𝑓(𝑡𝑛 , 𝑦𝑛 ) + 𝑓(𝑡𝑛+1 , 𝑦𝑛+1 )]
𝑡𝑛 2
𝐾1 𝐾2
Par introduction de la méthode d’Euler :
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓(𝑡𝑛 , 𝑦𝑛 ) 𝑜𝑢 ℎ 𝑒𝑠𝑡 𝑙𝑒 𝑝𝑎𝑠 𝑑′ 𝑖𝑛𝑡é𝑔𝑟𝑎𝑡𝑖𝑜𝑛

L’intégrale devient :
𝑘1 = 𝑓(𝑡𝑛 , 𝑦𝑛 )
ℎ
𝑦𝑛+1 − 𝑦𝑛 = (𝑘1 + 𝑘2 ) 𝑎𝑣𝑒𝑐 { 2 = 𝑓(𝑡𝑛 + ℎ, 𝑦𝑛 + ℎ𝑘1 )
𝑘
2
𝑦0 = 𝑦(𝑡0 )

D’où le schéma numérique pour l’intégration de l’équation différentielle :

ℎ
𝑦𝑛+1 = 𝑦𝑛 + (𝑘1 + 𝑘2 ) 𝑖 = 1 … 𝑛 𝑛: 𝑒𝑠𝑡 𝑙𝑒 𝑛𝑜𝑚𝑏𝑟𝑒 𝑑′𝑖𝑛𝑡𝑒𝑟𝑣𝑎𝑙𝑙𝑒𝑠
2

Algorithme
1- Définir h, t0, tmax
2- Définir f(t, y)
3- Définir le vecteur t allant de t0 à tmax par pas de h
4- Calcul de n nombre d’intervalles
5- Y(1)=y0
6- Boucle de compteur i sur les n intervalles
7- Calcul de k1 et k2
8- Calcul de y(i+1)
9- Fin de boucle
10- Tracer y(t)

La méthode de Runge-Kutta ordre 4

C’est une amélioration de méthode de Runge Kutta ordre 2. Elle la plus stable et la plus performante.

38
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Elle consiste estimer l’intégrale par la méthode de Simpson. La méthode de Simpson utilise un polynôme
d’interpolation sur trois points (voir chapitre sur intégration numérique). Pour un sous-intervalle i :
𝑡𝑛+1
ℎ
𝑦𝑛+1 − 𝑦𝑛 = ∫ 𝑓(𝑡, 𝑦)𝑑𝑡 = [𝑓(𝑡𝑛 , 𝑦𝑛 ) + 4𝑓(𝑡 1 , 𝑦 1 )]
𝑡𝑛 6 𝑛+
2
𝑛+
2

Le terme 4𝑓(𝑡𝑛+1 , 𝑦𝑛+1 ) peut-être écrit 2𝑓 (𝑡𝑛+1 , 𝑦𝑛+1 ) + 2𝑓(𝑡𝑛+1 , 𝑦𝑛+1 )

2 2 2 2 2 2

Introduisons ici la méthode d’Euler implicite pour le 1er terme et Euler explicite pour le 2eme terme :
h
𝑦 1 = 𝑦𝑖 + 𝑓(𝑡𝑖 , 𝑦𝑖 ) 𝑝𝑜𝑢𝑟 𝑙𝑒 1𝑒𝑟 𝑡𝑒𝑟𝑚𝑒
𝑖+
2 2
Et
h
𝑦 1 = 𝑦𝑖 + 𝑓 (𝑡 1 , 𝑦 1 ) 𝑝𝑜𝑢𝑟 𝑙𝑒 2𝑒𝑟 𝑡𝑒𝑟𝑚𝑒
𝑖+
2 2 𝑖+
2
𝑖+
2

Cela se traduit sur le schéma numérique suivant :

𝒌𝟏 = 𝒇(𝒕𝒏 , 𝒚𝒏 )
𝒉 𝒉
𝒌𝟐 = 𝒇 (𝒕𝒏 + , 𝒚𝒏 + 𝒌𝟏 )
𝒌𝟏 𝒌𝟐 𝒌𝟐 𝒌𝟒 𝟐 𝟐
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉 ( + + + ) 𝒂𝒗𝒆𝒄 𝒉 𝒉
𝟔 𝟑 𝟑 𝟔 𝒌𝟑 = 𝒇 (𝒕𝒏 + , 𝒚𝒏 + 𝒌𝟐 )
𝟐 𝟐
𝒌𝟒 = 𝒇(𝒕𝒏 + 𝒉, 𝒚𝒏 + 𝒉𝒌𝟑 )
{

Algorithme
1- Définir h, t0, tmax
2- Définir f(t, y)
3- Définir le vecteur t allant de t0 à tmax par pas de h
4- Calcul de n nombre d’intervalles
5- Y(1)=y0
6- Boucle de compteur i sur les n intervalles
7- Calcul de k1, k2, k3 et k4
8- Calcul de y(i+1)
9- Fin de boucle
10- Tracer y(t)

Travail demandé :

- Refaire les mêmes exercices que pour la méthode d’Euler

- Sur un même graphe, tracer les courbes de y(t) correspondant à chaque méthode
- Comparer les différentes méthodes avec la solution exacte. Et dire quelle est la méthode qui donne le
résultat le plus proche

39
USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

40
USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

PARTIE III
Petits projets d’application
Pour aller plus loin

41
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

La radiation du corps noir

Le corps noir est un corps matériel capable d’absorber tout le rayonnement électromagnétique qu’il reçoit
quel que soit la longueur d’onde. L’énergie du rayonnement absorbée aura pour effet d’agiter les constituants
du corps noir. Après la thermalisation, un rayonnement émis caractéristique est observé et un spectre d’émis-
sion qui ne dépend que de la température est ainsi observé.

Avant la célèbre loi de Planck de la distribution spectrale de la densité d’énergie du spectre d’émission du
corps noir, il y avait deux lois approximatives qui étaient répondus à l’époque à savoir :
 La loi de Raleigh-Jeans : valable pour les grandes longueurs d’onde.
 La loi de Wien : valable pour les petites longueurs d’onde.

La loi de Raleigh-Jeans :
8𝜋𝑘𝑇
𝑢(𝜆, 𝑇) =
𝜆4
La loi de Wien :

8𝜋ℎ𝑐 − ℎ𝑐
𝑢(𝜆, 𝑇) = 𝑒 𝜆𝑘𝑇
𝜆5
La loi de Planck :
8𝜋ℎ𝑐 1
𝑢(𝜆, 𝑇) = ℎ𝑐
𝜆5 −
𝑒 𝜆𝑘𝑇 −1

u : la densité d’énergie rayonnée par le corps noir pour la longueur d’onde λ. Son unité est : joule/m3

Dans cet exercice nous souhaitons faire un comparatif graphique de ces trois lois.

Pour cela on créera un nouveau script : corps_noir.m

On veut que à l’exécution du programme, l’utilisateur saisie la température à laquelle la comparaison sera
faite. Pour cela, on se sert de la fonction interne : input (‘chaine caractères’) qui met l’exécu-
tion en pause jusqu’à ce qu’une valeur soit saisie. Cette fonction retourne au programme la valeur de la tem-
pérature lue du clavier.

T = input ("Entrez la température en °K du corps noir : ") ;

Nous avons besoin de définir les constantes intervenantes dans les différentes formules (en unité MKSA mais
pas obligatoirement) :

c = 3*10^8; % la vitesse de la lumière dans le vide

h = 6.625*10.^-34; % la constante de Planck
k = 1.38*10.^-23; % la constante de Boltzmann

Pour discrétiser l’axe horizontal des longueurs d’onde ; nous allons utiliser un intervalle paramétré au lieu
d’un intervalle fixe pour remédier au problème du décalage du maximum de la courbe en fonction de la tem-
pérature. On choisira un intervalle [0, 5λmax] suffisamment large pour que la courbe de Raleigh-Jeans soit vi-
sible.

42
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Pour calculer la longueur d’onde correspondant au maximum ; on utilise par la formule de déplacement de
Wien : 𝜆𝑚𝑎𝑥 𝑇 = 0.0029 𝑚°𝐾 .

Lam0 = 0.0029/T ;
Lam = linspace(0,5*lam0,100); % 100 points dans L’intervalle de lambda
n = length(Lam); % n=100 nombre suffisant pour une courbe assez correcte

On définit la loi de Raileigh-Jeans :

Comme les opérations porte sur chaque élément du vecteur lambda on les précède par un . :

URJ= 8*pi*k*T ./ Lam.^4 ;

La loi de Wien :
UW = ((8*pi*h*c)./(Lam.^5)) .* exp(-h*c ./ (Lam*k*T));

La loi de Planck :
UP=(8*pi*h*c ./ Lam.^5) ./ (exp(h*c ./(k*T*Lam))-1);

La loi de Raileigh-Jeans devient très grande pour les petites longueurs d’onde les autres courbes ne seront
pas visibles devant elle. Pour remédier à cela, on doit limiter les courbes verticalement par à une valeur li-
mite. Cette valeur limite est prise à un de 1% au-dessus du maximum atteint par la courbe de Planck (qui est
presque la même que celle de Wien).

Umax = max(UP)*1.01; % limite verticale pour toutes les courbe

Tout est prêt pour réaliser les graphiques

plot(Lam,UP,'linewidth',2,'DisplayName','Planck''s law') ;
hold on
plot(Lam,UW,'--','linewidth',2,'DisplayName','Wien''s law');
plot(Lam,URJ,'.','linewidth',2,'DisplayName','Raileigh-Jeans''s law') ;
ylim([0 Umax*1.1]); % appliquer maintenant la limitation verticale

L’option ‘DisplayName’ permet de personnaliser le texte qui apparait dans la légende

legend('show');

grid on; % la grande grille

grid minor; % la grille fine

On affiche la température à droite du maximum de la courbe de Planck à une distance de 10% de λmax .
text(Lam0*1.1,Umax,strcat(num2str(T),' K'),'fontSize',14)

xlabel('\lambda(m)','fontsize',14)
ylabel('Energy density(J/m^4)','fontsize',14)
title('Blackbody Radiation','fontsize',14)

Pour afficher les limites du domaine visible ; on trace deux lignes à 0.4 et 0.8 μm avec les couleurs rouge et
bleue.

line([0.4 0.4]*1e-6, [0 Umax],'color','blue','linewidth',1);

line([0.8 0.8]*1e-6, [0 Umax],'color','red','linewidth',1);

43
 USTHB – Faculté de Physique - TP Analyse numérique - A. OUATIZERGA

Améliorations :

44