Matlab.pdf

Présentation de Matlab 1. Introduction - Historique 2. Démarrage de MATLAB 3. Génération de graphique avec MATLAB Systèmes d'équations linéaires 1. Matrices et Vecteurs dans MATLAB 2. Équations et systèmes linéaires dans MATLAB 3. Méthode directe (Méthode du pivot) 4....

Présentation de Matlab 1. Introduction - Historique 2. Démarrage de MATLAB 3. Génération de graphique avec MATLAB Systèmes d'équations linéaires 1. Matrices et Vecteurs dans MATLAB 2. Équations et systèmes linéaires dans MATLAB 3. Méthode directe (Méthode du pivot) 4. Méthodes itératives 4.1. Méthode de Jacobi 4.2. Méthode de Gauss-Seidel Polynômes et interpolation polynomiale Résolution des équations non linéaires 1. Opérations sur les polynômes dans MATLAB 1.1. Multiplication des polynômes 1.2. Division des polynômes 2. Manipulation de fonctions polynomiales dans MATLAB 2.1. Évaluation d’un polynôme 2.2. Interpolation au sens des moindres carrés 3. Interpolation linéaire et non linéaire 4. Interpolation de Lagrange 5. Résolution d’équations et de Systèmes d’équations non Linéaire 5.1. Résolution d’équations non Linéaires 5.2. Résolution de Systèmes d’équations non Linéaires Intégration numérique des fonctions 1. Introduction 2. Méthodes d’intégrations numériques 2.1. Méthode des trapèzes 2.2. Méthode de Simpson 3. Fonctions MATLAB utilisées pour l'intégration numérique Résolution numérique des équations différentielles et des équations aux dérivées partielles 1. Introduction 2. Équations différentielles du premier ordre 3. Équations différentielles du second ordre 4. Méthode de Runge-Kutta 4.1. Méthode de Runge-Kutta du second ordre 4.2. Méthode de Runge-Kutta à l'ordre 4 5. Méthode Matricielle avec des "Conditions aux Limites" 6. Conversion de coordonnées 6.1. Coordonnées polaires 6.2. Coordonnées cylindriques 6.3. Coordonnées sphériques 7. Problèmes en Coordonnées Cylindriques 8. Discrétisation de l'équation de la Conduction en régime instationnaire  Présentation de Matlab 1. Introduction - Historique MATLAB est une abréviation de Matrix LABoratory. Écrit à l’origine, en Fortran, par C. Moler, MATLAB était destiné à faciliter l’accès au logiciel matriciel développé dans les projets LINPACK et EISPACK. La version actuelle, écrite en C par the MathWorks Inc., existe en version professionnelle et en version étudiant. Sa disponibilité est assurée sur plusieurs plates- formes : Sun, Bull, HP, IBM, compatibles PC (DOS, Unix ou Windows), Macintoch, iMac et plusieurs machines parallèles. MATLAB est un environnement puissant, complet et facile à utiliser destiné au calcul scientifique. Il apporte aux ingénieurs, chercheurs et à tout scientifique un système interactif intégrant calcul numérique et visualisation. C'est un environnement performant, ouvert et programmable qui permet de remarquables gains de productivité et de créativité. MATLAB est un environnement complet, ouvert et extensible pour le calcul et la visualisation. Il dispose de plusieurs centaines (voire milliers, selon les versions et les modules optionnels autour du noyeau Matlab) de fonctions mathématiques, scientifiques et techniques. L'approche matricielle de MATLAB permet de traiter les données sans aucune limitation de taille et de réaliser des calculs numériques et symboliques de façon fiable et rapide. Grâce aux fonctions graphiques de MATLAB, il devient très facile de modifier interactivement les différents paramètres des graphiques pour les adapter selon nos souhaits. L'approche ouverte de MATLAB permet de construire un outil sur mesure. On peut inspecter le code source et les algorithmes des bibliothèques de fonctions (Toolboxes), modifier des fonctions existantes et ajouter d’autres. MATLAB possède son propre langage, intuitif et naturel qui permet des gains de temps de CPU spectaculaires par rapport à des langages comme le C, le TurboPascal et le Fortran. Avec MATLAB, on peut faire des liaisons de façon dynamique, à des programmes C ou Fortran, échanger des données avec d'autres applications (via la DDE : MATLAB serveur ou client) ou utiliser MATLAB comme moteur d'analyse et de visualisation. MATLAB comprend aussi un ensemble d'outils spécifiques à des domaines, appelés Toolboxes (ou Boîtes à Outils). Indispensables à la plupart des utilisateurs, les Boîtes à Outils sont des collections de fonctions qui étendent l'environnement MATLAB pour résoudre des catégories spécifiques de problèmes. Les domaines couverts sont très variés et comprennent notamment le traitement du signal, l'automatique, l'identification de systèmes, les réseaux de neurones, la logique floue, le calcul de structure, les statistiques, etc. MATLAB fait également partie d'un ensemble d'outils intégrés dédiés au Traitement du Signal. En complément du noyau de calcul MATLAB, l'environnement comprend des modules optionnels qui sont parfaitement intégrés à l'ensemble : 1) une vaste gamme de bibliothèques de fonctions spécialisées (Toolboxes) 2) Simulink, un environnement puissant de modélisation basée sur les schémas-blocs et de simulation de systèmes dynamiques linéaires et non linéaires 3) Des bibliothèques de blocs Simulink spécialisés (Blocksets) 4) D'autres modules dont un Compilateur, un générateur de code C, un accélérateur,... 5) Un ensemble d'outils intégrés dédiés au Traitement du Signal : le DSP Workshop. Quelles sont les particularités de MATLAB ? MATLAB permet le travail interactif soit en mode commande, soit en mode programmation ; tout en ayant toujours la possibilité de faire des visualisations graphiques. Considéré comme un des meilleurs langages de programmations (C ou Fortran), MATLAB possède les particularités suivantes par rapport à ces langages : la programmation facile, la continuité parmi les valeurs entières, réelles et complexes, la gamme étendue des nombres et leurs précisions, la bibliothèque mathématique très compréhensive, l’outil graphique qui inclus les fonctions d’interface graphique et les utilitaires, la possibilité de liaison avec les autres langages classiques de programmations (C ou Fortran). Dans MATLAB, aucune déclaration n’est à effectuer sur les nombres. En effet, il n'existe pas de distinction entre les nombres entiers, les nombres réels, les nombres complexes et la simple ou double précision. Cette caractéristique rend le mode de programmation très facile et très rapide. En Fortran par exemple, une subroutine est presque nécessaire pour chaque variable simple ou double précision, entière, réelle ou complexe. Dans MATLAB, aucune nécessité n’est demandée pour la séparation de ces variables. La bibliothèque des fonctions mathématiques dans MATLAB donne des analyses mathématiques très simples. En effet, l’utilisateur peut exécuter dans le mode commande n’importe quelle fonction mathématique se trouvant dans la bibliothèque sans avoir à recourir à la programmation. Pour l’interface graphique, des représentations scientifiques et même artistiques des objets peuvent être créées sur l’écran en utilisant les expressions mathématiques. Les graphiques sur MATLAB sont simples et attirent l’attention des utilisateurs, vu les possibilités importantes offertes par ce logiciel. MATLAB peut-il s’en passer de la nécessité de Fortran ou du C ? La réponse est non. En effet, le Fortran ou le C sont des langages importants pour les calculs de haute performance qui nécessitent une grande mémoire et un temps de calcul très long. Sans compilateur, les calculs sur MATLAB sont relativement lents par rapport au Fortran ou au C si les programmes comportent des boucles. Il est donc conseillé d'éviter les boucles, surtout si celles-ci est grande. 2. Démarrage de MATLAB Pour lancer l’exécution de MATLAB : sous Windows, il faut cliquer sur Démarrage, ensuite Programme, ensuite MATLAB, sous d’autres systèmes, se référer au manuel d’installation. L’invite ‘>>’ de MATLAB doit alors apparaître, à la suite duquel on entrera les commandes. La fonction "quit" permet de quitter MATLAB : >>quit La commande "help" permet de donner l’aide sur un problème donné. Exemple : >> help cos COS Cosine. COS(X) is the cosine of the elements of X. Autres commandes : what : liste les fichiers *.m et *.mat dans le directory utilisé who : liste les variables utilisées dans l’espace courant ans : réponse retournée après exécution d’une commande Exemple : >>x=[1:5,1] x= 123451 ou bien : >>[1:5,1] ans = 123451 clock : affiche l’année, le mois, le jour, l’heure, les minutes et les secondes. >>clock ans = 1.0e+003 * 1.9980 0.0100 0.0180 0.0170 0.0020 0.0098 >>date ans = 18-Oct-1998 Calcul en mode Commande dans MATLAB : Soit à calculer le volume suivant : où R=4cm Pour calculer V, on exécute les commandes suivantes : >>R=4 R= 4 >>V=4/3*pi*R^3 V= 268.0826 (Ici, pi=π). Calcul arithmétique : + plus - moins / division * multiplication Exemple : x=2, >>x=2 x= 2 >>P=(4*x^2-2*x+3)/(x^3+1) P= 1.6667 Test ‘if’ Ce test s'emploie, souvent, dans la plupart des programmes. Un test ‘if’ est toujours suivi par un ‘end’. Exemple : >>V=268.0826 V= 268.0826 >>if V>150, surface=pi*R^2, end surface = 50.2655 L’opérateur ‘NON’ Il est noté (ou symbolisé) par ‘~=’ Exemple : >>R=4 R= 4 >>if R~=2, V=4/3*pi*R^3;end L’opérateur ‘égal’ (==) dans ‘if’ Il est noté (ou symbolisé) par ‘==’. Exemple : >>R=4 R= 4 >>if R==4, V=4/3*pi*R^3;end L’opérateur ‘ou’ Il est noté (ou symbolisé) par ‘|’ Exemple : Si R=4 ou m=1, alors >>if R==4 | m==1, V=4/3*pi*R^3;end Autres opérateurs : > supérieur à < inférieur à >= supérieur ou égal supérieur à.* produit élément par élément de matrices < inférieur à.^ puissance élément par élément de matrices >= supérieur ou égal./ division élément par élément de matrices >A=zeros(n,m) où n et m sont respectivement le nombre de lignes et le nombre de colonnes. La matrice identité est définie dans Matlab par la fonction : ‘eye’ Exemple : C’est une matrice carrée d’ordre 5. Dans Matlab, on peut l’écrire : >>B=eye(5) B= 10000 01000 00100 00010 00001 Si Bt est la matrice transposée de B, cette dernière s’écrit dans Matlab : >>C=B' (l’appostrophe représente la transposée) C= 10000 01000 00100 00010 00001 Le déterminant d’une matrice carrée (C par exemple) est donné par : >>det(C) ans = 1 Soit A une matrice non nulle. La matrice inverse A-1 de A (si elle existe) est telle que : A*A-1=IId Dans Matlab, cette matrice inverse est donnée par : A-1=A^(-1)=inv(A)=pinv(A) Exemple : >>A=[1 3 5;2 -1 0;5 4 3] A= 135 2 -1 0 543 La matrice inverse de A est : >>inv(A) ans = -0.0682 0.2500 0.1136 -0.1364 -0.5000 0.2273 0.2955 0.2500 -0.1591 ou encore : >>pinv(A) ans = -0.0682 0.2500 0.1136 -0.1364 -0.5000 0.2273 0.2955 0.2500 -0.1591 2. Equations et systèmes linéaires dans MATLAB Considérons le systèmes suivants de n équations à m inconnues : Dans ce système, on a connue, connues et inconnues. Si n>m système sur-déterminé, Si n=m système de Cramer une solution unique (si det(A) ),  Si n>x=A^(-1)*y x= -0.5227 -1.0455 1.9318 Notation : Soit la matrice : Dans Matlab, la Matrice A peut être écrite sous les deux formes : >>A=[1,3,5;2,-1,0;5,4,3]; ou bien : >>A=[1 3 5;2 -1 0;5 4 3]; La séparation des lignes se fait par un point virgule et la séparation des colonnes se fait par un espace ou une virgule. Exemple : >>A=[1 3 5;2 -1 0;5 4 3] A= 135 2 -1 0 543 Problèmes mal conditionnés : Il existe un certain nombre d’équations linéaires solvables, mais leurs solutions sont incertaines à cause des erreurs d’arrondi lors des calculs. Les problèmes de ce type sont appelés ‘problèmes mal conditionnés’. L’arrondi durant les calculs ou les petits changements dans les coefficients calculés peuvent induire des erreurs significatives dans la résolution d’un problème (mal conditionné). Cet effet d’arrondi, peut être illustré par le système des 2 équations suivantes, par exemple : La solution de ce système est : Pour monter l’erreur introduite sur les solutions en modifiant les coefficients (arrondi), on donne au coefficient ‘2,01045’ de la première équation une légère augmentation de ‘0,001’. Dans ce cas, le système devient : Ceci met en évidence l'erreur engendrée par un arrondi. La matrice C correspondant au système précédent est : La norme de C s’écrit : Dans Matlab on définit le nombre-condition de la matrice C [notation cond(C)] par : Ce nombre satisfait toujours la condition : Pour les matrices inconditionnées, ce nombre-condition attiend des grandes valeurs, mais ne donne pas une estimation directe de l’erreur sur la solution. Dans Matlab, ce nombre se calcule par la fonction ‘cond(C)’ où C est la matrice. Exemple1 : >>C C= 0.1206 0.9878 0.9876 >>cond(C) ans = 6.5598e+003 Exemple2 : L'exemple de la matrice d’Hilbert représente un problème à matrice mal conditionnée. Cette matrice est définie par : , où Programmer le nombre-condition (cond(A)) ainsi que det(A).det(A-1) pour une matrice d’Hilbert de 5 par 5 à 14 par 14 : Solution : clear; for n=5:14 for i=1:n for j=1:n a(i,j)=1/(i+j-1); end end C=cond(a); d=det(a)*det(a^(-1)); fprintf('n=%3.0f\t cond(a)=%e\t det*det=%e\n',n,C,d); end Après exécution du programme, on obtient les résultats suivants : n= 5 cond(a)=4.766073e+005 det*det=1.000000e+000 n= 6 cond(a)=1.495106e+007 det*det=1.000000e+000 n= 7 cond(a)=4.753674e+008 det*det=1.000000e+000 n= 8 cond(a)=1.525758e+010 det*det=9.999999e-001 n= 9 cond(a)=4.931532e+011 det*det=1.000001e+000 n= 10 cond(a)=1.602534e+013 det*det=9.999812e-001 n= 11 cond(a)=5.218389e+014 det*det=1.000119e+000 n= 12 cond(a)=1.768065e+016 det*det=1.015201e+000 n= 13 cond(a)=3.682278e+018 det*det=9.707419e-001 n= 14 cond(a)=1.557018e+018 det*det=-4.325761e+000 Durant l’exécution du programme, un message s’affiche sur l’écran pour n=11 à n=14 : Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.659249e-017. Ce qui signifie qu’à partir de n=11 jusqu’à n=14, on a det(a)*det(a-1) différent de 1, donc les erreurs d’arrondi commencent leurs apparitions. 3. Méthode directe (Méthode du pivot) Soit à résoudre le système suivant de 3 équations à 3 inconnues par la méthode du pivot (dans Matlab) : déterminer. Solution : 1. On définit tout d’abord la matrice argument dans Matlab par : >>A=[-0.04 0.04 0.12 3;0.56 -1.56 0.32 1; -0.24 1.24 -0.28 0] ; 2. On choisit la ligne 1 comme ligne pivot : a(1,1)=-0.04, on divise la ligne 1 par a(1,1) b(1,:)=a(1,:)/a(1,1) nouvelle ligne, on annule le 1er terme de la ligne 2 et le 1er terme de la ligne 3 : b(2,:)=a(2,:)-b(1,:)*a(2,1), b(3,:)=a(3,:)-b(1,:)*a(3,1). On obtient après calculs : b= 1 -1 -3 -75 0 -1 2 43 0 1 -1 -18 3. On choisit la ligne 2 comme ligne pivot. On divise cette ligne par b(2,2)=-1, on obtient : c(2,:)=b(2,:)/b(2,2) nouvelle ligne = 0 1 -2 -43 on annule le terme b(1,2) de la ligne 1 et le terme b(3.2) de la ligne 3 : c(1,:)=b(1,:)-c(2,:)*b(1,2) 1 0 -5 118, c(3,:)=b(3,:)-c(2,:)*b(3,2) 0 0 1 25. 4. On choisit la ligne 3 comme ligne pivot c(3, :). On divise cette ligne par c(3,3)=1, on obtient : d(3,:)=c(3,:) nouvelle ligne = 0 0 1 25 on annule dans la ligne 1 c(1,3) et dans la ligne 2 c(2,3) :  d(1,:)=c(1,:)-d(3,:)*c(1,3) 1 0 0 7, d(2,:)=c(2,:)-d(3,:)*c(2,3) 0 1 0 7. D’où la matrice d s’écrit : d= 1007 0107 0 0 1 25 et la solution est : L’algorithme suivant (dans Matlab) permet de résoudre le système précédent : clear %effacer toutes les variables en mémoire dans Matlab a=[-0.04 0.04 0.12 3;0.56 -1.56 0.32 1;-0.24 1.24 -0.28 0];a x=[0 0 0]'; %x est le vecteur colonne de composantes xi qui est initialisé ici % 1er point pivot ligne 1 -calculs sur les lignes 2 et 3 b(1,:)=a(1,:)/a(1,1); b(2,:)=a(2,:)-b(1,:)*a(2,1); b(3,:)=a(3,:)-b(1,:)*a(3,1);b % 2ème pivot ligne 2 - calculs sur les lignes 1 et 3 c(2,:)=b(2,:)/b(2,2); c(1,:)=b(1,:)-c(2,:)*b(1,2); c(3,:)=b(3,:)-c(2,:)*b(3,2);c % 3ème pivot ligne 3 - calculs sur les lignes 1 et 2 d(3,:)=c(3,:)/c(3,3); d(1,:)=c(1,:)-d(3,:)*c(1,3); d(2,:)=c(2,:)-d(3,:)*c(2,3);d %Solutions recherchées x(1)=d(1,4); x(2)=d(2,4); x(3)=d(3,4);x après l’exécution de ce programme, on obtient : a= -0.0400 0.0400 0.1200 3.0000 0.5600 -1.5600 0.3200 1.0000 -0.2400 1.2400 -0.2800 0 b= 1.0000 -1.0000 -3.0000 -75.0000 0 -1.0000 2.0000 43.0000 0 1.0000 -1.0000 -18.0000 c= 1.0000 0 -5.0000 -118.0000 0 1.0000 -2.0000 -43.0000 0 0 1.0000 25.0000 d= 1.0000 0 0 7.0000 0 1.0000 0 7.0000 0 0 1.0000 25.0000 x= 7.0000 7.0000 25.0000 4. Méthodes itératives 4.1. Méthode de Jacobi Résoudre le système suivant par la méthode itérative de Jacobi : Écrire un algorithme permettant de calculer x1, x2 et x3 par itérations. Déterminer le nombre d’itérations nécessaires pour obtenir :. Solution : Appelons l’algorithme permettant de calculer la solution ‘jacobi.m’. Dans cet algorithme, nous avons utilisé un facteur de relaxation "λ" pour réaliser la convergence, car sinon le système diverge. 4.2. Méthode de Gauss-Seidel Résoudre le même système précédent par la méthode de Gauss-Seidel. Écrire un algorithme (Matlab) permettant de résoudre ce système et déterminer le nombre d’itérations nécessaires pour obtenir une erreur. Comparer les résultats avec la méthode de Jacobi et conclure. Représenter sur le même graphe l’évolution des solutions en fonction du nombre d’itérations. Changer les valeurs initiales , exécuter le programme puis conclure. Solution : * Méthode itérative générale : Appelons l’algorithme permettant de calculer la solution ‘G-seid.m’ On se donne des valeurs arbitraires : par exemple =0. Le système à résoudre est : A*x=y ; où , et La matrice A et le vecteur colonne y sont connus, le vecteur colonne x reste inconnu. Ainsi, la méthode itérative générale consiste à écrire : pour et. En particulier, pour un système de 3 équations à 3 inconnues, on a : i=1 j=2 et 3 i=2 j=1 et 3  i=3 j=1 et 2 et ceci pour chaque itération. Programme itmg.m : %****************************** % Résolution d'un système linéaire * % par la méthode itérative générale * %****************************** clear all; clc; fprintf('Méthode itérative générale\n'); n=30000; a=[1 1 1;2 -1 3;3 2 -2]; y=[1 4 -2]; x=zeros(1,3); w=0.2; % facteur de relaxation : 0>f=[3 2 -1 4]; >>g=[2 0 -3 5 -1]; >>h=conv(f,g) h= 6 4 -11 17 10 -19 21 -4 Ainsi, le polynôme obtenu est : 1.2. Division des polynômes La fonction ‘deconv’ donne le rapport de convolution de deux polynômes (déconvolution des coefficients du polynôme). L’exemple suivant montre l’utilisation de cette fonction. Soient les mêmes fonctions précédentes et : La division de par : est donnée par la fonction ‘deconv’ : >>f=[3 2 -1 4]; >>g=[2 0 -3 5 -1]; >>h=deconv(g,f) h= 0.6667 -0.4444 et le polynôme obtenu est : 2. Manipulation de fonctions polynomiales dans MATLAB Soit le polynôme suivant : où n est degré du polynôme et sont les coefficients du polynôme. Ce polynôme peut être écrit sous la forme : Après factorisation, on a : où sont les racines du polynômes. Exemple : Ce polynôme est équivalent à : ou encore : Un polynôme d’ordre n possède n racines qui peuvent être réelles ou complexes. Dans MATLAB, un polynôme est représenté par un vecteur contenant les coefficients dans un ordre décroissant. Exemple : Le polynôme : qui est représenté dans MATLAB par : >>P=[2 1 4 5]; a pour racines. Pour trouver ces racines, on doit exécuter la fonction ‘roots’. D’où : >>r=roots(P); et le résultat donné est : >> r r= 0.2500 + 1.5612i 0.2500 - 1.5612i -1.0000 Les trois racines de ce polynôme (dont 2 sont complexes) sont données sous forme d’un vecteur colonne. Quand les racines sont connues, les coefficients peuvent être recalculés par la commande ‘poly’. Exemple : >>poly(r) ans = 1.0000 0.5000 2.0000 2.5000 La fonction ‘poly’ accepte aussi une matrice comme argument dont elle retourne le polynôme caractéristique. Exemple : >>A=[3 1;2 4]; >>p=poly(A) p= 1 -7 10 Ainsi, le polynôme caractéristique de la matrice A est : Les racines de ce polynôme sont les valeurs propres de la matrice A. ces racines peuvent être obtenues par la fonction ‘eig’ : >>Val_prop=eig(A) Val_prop = 2 5 2.1. Évaluation d’un polynôme Pour évaluer le polynôme en un point donné, on doit utiliser la fonction ‘polyval’. On évalue ce polynôme pour x=1, par exemple : >>polyval(P,1) ans = 12 Exemple : On veut évaluer en x=2.5 le polynôme suivant : >>C=[3 -7 2 1 1]; >>x=2.5; >> y=polyval(C,x) y= 23.8125 Si x est un vecteur contenant plusieurs valeurs, y sera aussi un vecteur qui contiendra le même nombre d’éléments que x. 2.2. Interpolation au sens des moindres carrés Dans le domaine de l’analyse numérique des données, on a souvent besoin d’établir un modèle mathématique liant plusieurs séries de données expérimentales. L’interpolation polynomiale consiste à approcher la courbe liant les deux séries de mesures par un polynôme. Les coefficients optimaux de ce polynôme sont ceux qui minimisent la variance de l’erreur d’interpolation. Ce principe (voir cours) est connu sous le nom de la méthode des moindres carrés. La fonction ‘polyfit’ retourne le polynôme P de degré n permettant d’approcher la courbe au sens des moindres carrés. Exemple : >> x=[1.1 2.3 3.9 5.1]; >>y=[3.887 4.276 4.651 2.117]; >>a=polyfit(x,y,length(x)-1) a= -0.2015 1.4385 -2.7477 5.4370 Ainsi, le polynôme d’interpolation de y (d’ordre length(x)-1=3) est : Pour déduire l’erreur entre les valeurs expérimentales et le modèle obtenu par la fonction ‘polyfit’, on dispose de la fonction ‘polyval’ qui retourne la valeur du polynôme P pour toutes les composantes du vecteur (ou de la matrice) x. Ainsi, cette fonction donne : >>yi=polyval(a,x) yi = 3.8870 4.2760 4.6510 2.1170 On remarque ici que les valeurs expérimentales de y sont bien restituées (y=yi). Dans ce cas, on a un coefficient de corrélation qui est égal à 1 (voir la définition de ce coefficient dans le cours). Pour mieux comprendre ces fonctions prédéfinis dans MATLAB, on va simuler une courbe expérimentale par une sigmoïde à laquelle on superpose un bruit du type Gaussien. Cette courbe sera donnée par : Le programme correspondant à l’approximation des ces données dans MATLAB (regres.m, par exemple) est le suivant : %****************************************** % Génération de données expérimentales: * % y=1./(1+exp(-x))+0.05*randn(1,length(x))* %****************************************** clc; % Effacer l'écran clear all; % Effacer des variables de l'espace de travail x=-5:0.1:5; % Intervalle de définition et de calcul de la sigmoïde % Fonction sigmoïde bruitée y=1./(1+exp(-x))+0.05*randn(1,length(x)); plot (x,y); % Tracé de la sigmoïde bruitée title('Fonction sigmoïde bruitée - Polynôme d''interpolation'); xlabel('x');ylabel('y'); % Polynôme d'interpolation d'ordre 1 P=polyfit(x,y,1); % Valeurs du polynôme d'interpolation Vp=polyval(P,x); % Tracé du polynôme d'interpolation hold on; plot(x,Vp,'--'); % Calcul de l'erreur d'interpolation erreur=y-Vp; % Tracé de la courbe de l'erreur plot(x,erreur,':') grid gtext('Mesures') gtext('Erreur') gtext('Modèle') hold off % Affichage du polynôme d'interpolation disp('Polynôme d''interpolation') P Var_erreur=num2str(std(erreur).^2); disp(['La variance de l''erreur d''interpolation est : ',Var_erreur]) Après exécution du programme, on obtient à l’ordre 1 (droite affine : fig. 1 ci-dessous) les résultats suivants : >> regres Polynôme d'interpolation P= 0.1309 0.5008 La variance de l'erreur d'interpolation est : 0.011277 On remarque ici que le polynôme d’interpolation d’ordre 1 n’est pas une bonne approximation pour ces données. En effet, un polynôme d’ordre 5 donne une meilleure approximation de la sigmoïde (fig. 2). Ainsi, dans le programme précédent (sigreg.m), on change dans le ‘1’ par ‘5’ dans la fonction ‘polyfit’, et on obtient les résultats suivants : >> regres Polynôme d'interpolation P= 0.0002 -0.0000 -0.0111 0.0008 0.2326 0.4844 La variance de l'erreur d'interpolation est : 0.002279 Fig.1 : Interpolation linéaire  Fig.2 : Interpolation par un polynôme d’ordre 5 3. Interpolation linéaire et non linéaire Une interpolation consiste à relier les points expérimentaux par une courbe sous forme de segments de droites ou de courbes polynomiales. Ceci peut être réalisé par la fonction ‘interp1’. La commande ‘interp1(x,y,xi,’type’)’ retourne un vecteur de mêmes dimensions que xi et dont les valeurs correspondent aux images des éléments de xi déterminées par interpolation sur x et y. Si f est l’interpolation de y, la chaîne ‘type’ spécifie alors le type d’interpolation qui doit être parmi les suivants : ‘linear’ interpolation linéaire ‘spline’ interpolation par splines cubiques, ‘cubic’ interpolation cubique. Si on ne spécifie pas le type, l’interpolation linéaire est choisie par défaut. Exemple 1 : >>x = 0:10; y = sin(x); xi = 0:.25:10; >>yi = interp1(x,y,xi,’cubic’); plot(x,y,'o',xi,yi)  Fig. 3 : Interpolation cubique Exemple 2 : Dans le cas suivant, on étudiera ces différents types d’interpolation sur un même exemple de valeurs discrètes de la fonction ‘cosinus’. On appellera l’algorithme : ‘interpol.m’. %************************************* % Utilisation de la fonction interp1 * %************************************* clear all; clc; x=0:10; y=cos(x); % Points à interpoler z=0:0.25:10; % Le pas du vecteur z est inférieur à celui de x % Interpolation linéaire figure(1); f=interp1(x,y,z); % Tracé des valeurs réelles et de la courbe d'interpolation plot(x,y,'*r',z,f); grid on; xlabel('Interpolation'); % Interpolation par splines cubiques figure(2); f=interp1(x,y,z,'spline'); plot(x,y,'*r',z,f); grid on; xlabel('Interpolation par splines cubiques'); En exécutant ce programme, on obtient les courbes suivantes : Fig. 4 : Interpolation linéaire  Fig. 5 : Interpolation par splines cubiques La fonction ‘interp2’ réalise l’interpolation dans l’espace trois dimensions (3D). 4. Interpolation de Lagrange Exemple : Les masses volumiques du matériau pour différentes températures sont données par le tableau ci- dessous : i 1 2 3 Température T (en °C) 94 205 371 Masse volumique R(T) : (en kg/m3) 929 902 860 Écrire la formule d’interpolation de Lagrange qui permet d’interpoler les différents points de données précédentes. Trouver les masses volumiques du sodium pour T=251 °C, T=305 °C et T=800 °C en utilisant l'interpolation de Lagrange. Solution : a) Puisque le nombre de points est égal à 3, alors le polynôme de Lagrange sera de degré 2. Ce polynôme s’écrit : Soit : D’où pour T=251 °C, on obtient R(251) = 890.5 kg/m3. L'algorithme de calcul de R(251), R(305) et R(800) est : ‘polylag.m’ qui est listé ci-dessous : %**************************************** % Interpolation polynomiale de Lagrange * %**************************************** clc; clear; T=[94 205 371]; R=[929 902 860]; Ti=[251 305 800]; Ri=lag(T,R,Ti) La fonction ‘lag.m’ est un sous programme permettant de donner l’interpolation de Lagrange. La liste de ce sous programme est : %*************************** % Sous-programme Lagran_.m * %*************************** function Ri=lag(T,R,Ti) Ri=zeros(size(Ti)); % Initialisation des Ti n=length(R); % Nombre de points for i=1:n y=ones(size(Ti)); for j=1:n if i~=j y=y.*(Ti-T(j))/(T(i)-T(j)); end Ri=Ri+y*R(i) end end return Il suffit maintenant d’exécuter le programme ‘polylag.m’ pour obtenir les résultats de Ri. Ces résultats sont les suivants : >>Ri Ri = 890.5561 876.9316 742.45559 5. Résolution d’équations et de Systèmes d’équations non Linéaire 5.1. Résolution d’équations non Linéaires Dans cette partie, on s’intéressera à la méthode de Newton-Raphson pour la résolution des équations non linéaires à une ou à plusieurs variables. On cherchera ensuite à trouver la valeur qui annulera la fonction. La méthode de Newton-Raphson permet d’approcher par itérations la valeur au moyen de la relation suivante : Dans toutes les méthodes itératives, il est nécessaire pour éviter une divergence de la solution, de bien choisir la valeur initiale. Celle-ci peut être obtenue graphiquement. On se propose d’appliquer cette méthode pour la recherche des racines de la fonction non linéaire suivante : Dans un premier temps, on se propose de tracer la courbe représentative de cette fonction en utilisant le programme ci-dessous ‘pnum1.m’: %************************* % Etude de la fonction : * % f(x)=exp(x)-2.cos(x) * %************************* x=-1:0.1:1; f=exp(x)-2*cos(x); plot(x,f); grid on; title('Fonction : f(x)=exp(x)-2.cos(x)'); Après exécution du programme, on obtient la courbe sur la figure 6. D ‘après cette courbe, il est judicieux de choisir un ; car est proche de zéro pour avoir une convergence rapide. La fonction dérivée a pour expression :. fig. 6 : Localisation de la valeur où Pour chercher la solution de , on peut rajouter au programme précédent ‘pnum1.m’ quelques lignes : %************************* % Etude de la fonction : * % f(x)=exp(x)-2.cos(x) * %************************* clf; x=-1:0.1:1; f=exp(x)-2*cos(x); figure(1); plot(x,f); grid on; title('Fonction : f(x)=exp(x)-2.cos(x)'); clear all; clc; x(1)=input('Donner la valeur initiale x(1): \n'); e=1e-10; n=5000; for i=2:n f=exp(x(i-1))-2*cos(x(i-1)); diff=exp(x(i-1))+2*sin(x(i-1)); x(i)=x(i-1)-f/diff; if abs(x(i)-x(i-1))>x x= 0.8500 0.7800 0.7719 0.7718 0.7718 0.7718 0.7718 0.7718 0.7718 0.7718 >>y y= 0.3500 0.4271 0.4197 0.4196 0.4196 0.4196 0.4196 0.4196 0.4196 0.4196 Les 2 solutions peuvent être obtenues au bout de 2 itérations grâce au bon choix des conditions initiales. Avec les conditions initiales x(1)=0,2 et y(1)=0.1, il y a convergence vers la solution (0,0) au bout de 3 itérations.  Fig. 11 : Convergence des solutions (x,y) : x(1)=0,85 et y(1)=0,35 Pour les conditions initiales x(1)=0,2 et y(1)=0,1, l’évolution des solutions x et y est : >>x x= 0.2000 -0.1031 -0.0112 -0.0002 -0.0000 -0.0000 -0.0000 0 0 0 >>y y= 0.1000 -0.0594 -0.0054 -0.0001 -0.0000 -0.0000 -0.0000 0 0 0 La boite à outil ‘Optimisation Toolbox’ propose les fonctions ‘fsolve’ et ‘fsolve2’ qui permettent respectivement, la recherche des racines d’une fonction dans un intervalle donné et les solutions d’équations non linéaires et de systèmes d’équations non linéaires. Fig. 12 : Convergence des solutions (x,y) : x(1)=0,2 et y(1)=0,1  Intégration numérique des fonctions 1. Introduction Dans la plupart des cas, les fonctions analytiques, du fait de leurs complexités, ne sont pas intégrables analytiquement. Dans d'autres cas, on a des fonctions qui sont évaluées numériquement en différents points de l’intervalle où ces dernières sont données, et l’intégrale de ces types de fonctions ne peut être obtenue que par des approches numériques. Dans ce chapitre, on s’intéresse aux méthodes utilisées fréquemment ; à savoir la méthode des trapèzes, la méthode de Simpson, les formules de Newton-Cotes et la méthode de Gauss. Nous étudierons également les méthodes de calcul d’intégrales doubles et les intégrales des fonctions non définies sur leurs bornes. 2. Méthodes d’intégrations numériques 2.1. Méthode des trapèzes Soit la fonction à intégrer sur. L’intégrale de s’écrit en utilisant la méthode des trapèzes : où ; ; et. Le terme représentant l’erreur est : est la moyenne de sur l’intervalle. L’erreur est inversement proportionnelle à la valeur de. Si la fonction est donnée sur les intervalles réguliers ( ), la méthode des trapèzes peut être écrite dans Matlab sous la forme : avec : On peut également faire un programme pour calculer l'intégrale. Celui-ci appelé 'trapez_v.m' par exemple, est listé ci-dessous : function I=trapez_v(g,h) I=(sum(f)-(f(1)+f(length(f)))/2)*h; Considérons par exemple la fonction à intégrer : sur un intervalle où le pas est "h" égal à 1. En mode interactif dans MATLAB (mode commande), avant de lancer le programme 'trapez_v.m', on donne la fonction ainsi que ses bornes : >>x=-10:1:8; >>f=x.^2+2*x-1; >>h=1; On exécute ensuite le programme 'trapez_v.m', et on obtient : >>I=trapez_v(f,h) I= 2265 En utilisant les programmes 'trapez_n.m' et 'trapez_g.m' listés ci-après, on peut calculer aussi l'intégrale I. >>I=trapez_n('fonction_f',a,b,n) ou >>I=trapez_g('fonction_f',a,b,n) Liste du programme 'trapez_n.m' : function I=trapez_n(fonction_f,a,b,n) h=(b-a)/n; x=a+(0:n)*h; f=feval(fonction_f,x); I=trapez_v(f,h) Liste du programme 'trapez_g.m' : function I=trapez_g(fonction_f,a,b,n); n=n; hold off; h=(b-a)/n; x=a+(0:n)*h; f=feval(‘fonction_f’,x); I=h/2*(f(1)+f(n+1)); if n>1 I=I+sum(f(2:n))*h; end h2=(b-a)/100; xc=a+(0:100)*h2; fc=feval(‘fonction_f’,xc); plot(xc,fc,'r'); hold on; title('Méthode des trapèzes'); xlabel('x'); ylabel('y'); grid on; plot(x,f,'m'); plot(x,zeros(size(x)),'c') for i=1:n; plot([x(i),x(i)],[0,f(i)],'g'); end La fonction 'fonction_f' est donnée par le programme 'fonction_f.m'. Remarque : La fonction y=feval('sin',x) est équivalente à y=sin(x). Exemples : 1°) On donne la fonction 'fonction_f' par le sous-programme 'fonction_f.m' listé ci-dessous : %**************** % fonction f(x) * %**************** function f=f(x); f=x+2*log(0.0000025235*x); En exécutant le programme trapez_g('fonction_f',1,50,60), où les bornes et sont égales respectivement à 1 et 50 et le nombre d'intervalles est égal 60, on obtient l'intégrale de sur pour 60 pas : >>trapez_g('fonction_f',1,50,60) ans = 279.3889  Fig. 1 : Intégrale de la fonction 2°) Un véhicule de masse m=2000 kg se déplace à la vitesse de V=30 m/s. Soudainement, le moteur est débrayé à. L'équation du mouvement après l'instant est donnée par : (a) où x est la distance linéaire mesurée à partir de. Dans cette équation (a), le terme de gauche représente la force d'accélération, le 1er terme de droite représente la résistance aérodynamique exercée par le vent sur le véhicule et le second terme est le coefficient de frottement. Calculer la distance au delà de laquelle la vitesse de la voiture se réduit à 15 m/s. Solution : L'équation (a) peut être écrite sous la forme : L'intégration de cette équation donne : Cette intégrale peut être évaluée par la méthode des trapèzes. Si on se donne 60 intervalles (ou 61 points), on peut écrire : où et. En définissant : et en appliquant la méthode des trapèzes pour l'intégration, on obtient : (b) Le programme suivant (trapeze.m)permet de calculer la valeur donnée par la formulation précédente (b) : >>trapeze x= 127.5066 En comparant cette solution à la solution exacte (=127,51 m), on remarque que l'erreur relative induite par cette méthode est de l'ordre de 2,7.10-3 %. La liste du programme trapeze.m est la suivante : clear; clc; nb_points=61; i=1:nb_points; h=(30-15)/(nb_points-1); V=15+(i-1)*h; f=2000*V./(8.1*V.^2+1200); x=trapez_v(f,h) 3°) Connaissant l'intégrale exacte de la fonction sur l'intervalle , on a. Étudier l'effet du nombre d'intervalles sur l'erreur de calcul en adoptant la méthode trapézoïdale : Solution : Le programme suivant 'trz.m' permet de calculer l'intégrale : clc; clear; Iexact=4.006994; a=0;b=2; fprintf('\n Méthode Trapézoïdale étendue \n'); fprintf('\n n\t \t I\t Pourc. erreur relat.\n'); fprintf('---------------------------------------------\n'); n=1; for k=1:10 n=2*n; h=(b-a)/n; i=1:n+1; x=a+(i-1)*h; f=sqrt(1+exp(x)); I=trapez_v(f,h); erreur=abs(Iexact-I)/Iexact; fprintf('%d\t %10.5f\t %10.8f\n',n,I,erreur); end Les résultats obtenus sont les suivants : Méthode Trapézoïdale étendue n I Pourc. erreur relat. --------------------------------------------- 2 4.08358 0.01911429 4 4.02619 0.00478998 8 4.01180 0.00119825 16 4.00819 0.00029965 32 4.00729 0.00007496 64 4.00707 0.00001878 128 4.00701 0.00000474 256 4.00700 0.00000123 512 4.00700 0.00000035 1024 4.00699 0.00000013 D'après ces résultats, on montre que quand est doublé, l'erreur relative décroît par un facteur 4. 2.2. Méthode de Simpson Soit l'intégrale de sur l'intervalle. Par la méthode de Simpson, s'écrit : où ; et. Le terme est donné par : et est la moyenne de sur l'intervalle ouvert. Exemple : Évaluer l'intégrale de sur l'intervalle avec la méthode de Simpson pour et. Solution : La liste du programme 'msimp.m' est la suivante : clc; clear; Iexact=4.006994; a=0;b=2; fprintf('\n Méthode de Simpson \n'); fprintf('\n n\t \t I\t Pourc. erreur relat.\n'); fprintf('---------------------------------------------\n'); n=1; for k=1:4 n=2*n; h=(b-a)/n; i=1:n+1; x=a+(i-1)*h; f=sqrt(1+exp(x)); I=h/3*(f(1)+4*sum(f(2:2:n))+f(n+1)); if n>2 I=I+h/3*2*sum(f(3:2:n)); end  erreur=abs(Iexact-I)/Iexact; fprintf('%d\t %10.5f\t %10.8f\n',n,I,erreur); end Les résultats de cet algorithme sont donnés ci-dessous : Méthode de Simpson n I Pourc. erreur relat. --------------------------------------------- 2 4.00791 0.00022935 4 4.00705 0.00001521 8 4.00700 0.00000101 16 4.00699 0.00000012 On montre qu'à partir de , l'erreur relative devient presque nulle (de l'ordre de ). Les programmes 'simp_v' et 'simp_n' peuvent être utilisés pour intégrer la fonction. Dans l'algorithme 'simp_v(f,h)', est un vecteur contenant les ordonnées de l'intégrante et représente le pas de l'intervalle d'intégration. D'où : >>I=simp_n('fonction_f',a,b,n) permet d'évaluer l'intégrale dans laquelle 'fonction_f' est le nom du sous-programme contenant la fonction , et sont les limites de l'intervalle d'intégration et. Exemple : Soit où et. Ecrire un algorithme utilisant le sous-programme simp_v(f,h) (à écrire) et qui permet d'évaluer sur 100 intervalles. Solution : Appelons ce programme simp.m par exemple, et donnons sa liste : clear; R=5;g=9.81;b=0.1; x1=-0.90*R;x2=R; h=(x2-x1)/100; x=x1:h:x2; f=(R^2-x.^2)./(b^2*sqrt(2*g*(x+R))); I=simp_v(f,h) Le sous-programme simp_v(f,h) est : function I=simp_v(f,h); n=length(f)-1; % Nombre d’intervalles if n==1 fprintf('Données à un seul intervalle'\n'); return; end if n==2 I=h/3*(f(1)+4*f(2)+f(3)); return; % Revenir et continuer le calcul end if n==3 I=(3/8)*h*(f(1)+3*f(2)+3*f(3)+f(4)); return; end I=0; if 2*floor(n/2)~=n ; %Le nombre d’intervalle doit être pair I=3*h/8*(f(n-2)+3*f(n-1)+3*f(n)+f(n+1)); m=n-3; else m=n; end I=I+h/3*(f(1)+4*sum(f(2:2:m))+f(m+1)); if m>2 I=I+(h/3)*2*sum(f(3:2:m)) end En exécutant le programme simp.m, on trouve l'intégrale par la méthode de Simpson.: >>simp I= 1.8522e+003 3. Fonctions MATLAB utilisées pour l'intégration numérique Dans Matlab (Toolbox), il existe 2 fonctions appelées 'quad' et 'quad8' pour l'intégration numérique. La fonction 'quad' utilise la méthode de Simpson et la fonction 'quad8' utilise les formules de Newton-Cotes à l'ordre 8. Ces deux fonctions sont utilisées de la façon suivante : quad('fonction_f',a,b) quad('fonction_f',a,b,tol) quad('fonction_f',a,b,tol,trace) Dans la première forme, la tolérance 'tol' qui correspond à l'erreur relative , est considérée égale à 0,001 par défaut. Ainsi, le calcul quadratique de l'intégrale est réitéré jusqu'à ce que la tolérance soit satisfaite. Si la 3ème forme est utilisée avec une valeur non nulle de 'trace', un graphique d'évolution des itérations sera affiché sur l'écran. Quand à la fonction 'quad8', elle est utilisée de la même manière que la fonction 'quad' dans les 3 formes d'expressions précédentes. Exemple : Soit à calculer l'intégrale suivante : avec la transformation suivante : Dans ce cas, on a : On obtient donc : Or, on connaît la fonction (fonction Gamma) eulérienne du deuxième espèce qui est définie par : où est un réel. Pour un entier strictement positif, a la propriété suivante : On pourra donc s'en servir pour calculer la factorielle d'un entier naturel : Dans Matlab, cette fonction est notée 'gamma'. Exemple : La factorielle de 10 est donnée par gamma(10+1)=gamma(11) : >> fact=gamma(11) fact = 3628800 Cette propriété est à l'origine du nom de la 'fonction factorielle', attribuée souvent à la fonction 'Gamma'. La fonction 'gamma' est donc prédéfinie dans Matlab. Cette fonction permet, de part sa définition pour tous les arguments réels positifs, d'introduire 'la factorielle' des arguments non entiers. Exemple : >> format long >>gamma(0.5) ans = 1.77245385090552 >>racine_de_pi=sqrt(pi) racine_de_pi = 1.77245385090552 >>gamma(0) Warning: Divide by zero. ans = Inf >>gamma(-2) Warning: Divide by zero. ans = -Inf Tracé de la fonction Gamma Tracer la fonction Gamma pour. >>x=-3:0.1:3; >>gx=gamma(x); >>plot(x,gx);grid;title('Allure de la fonction Gamma'); >>xlabel('x');ylabel('Gamma(x)'); >>hold on;plot(x,zeros(size(x))) Fig. 2 : Fonction Revenons maintenant à l'intégration de la fonction : D'après le changement de variable effectué, cette intégrale s'écrit : Calculons cette intégrale par les fonctions Matlab 'quad8' et 'gamma'. La fonction 'quad8' nécessite l'écriture de la fonction à intégrer dans un fichier contenant l'extension '.m' (programme Matlab). Liste du programme fct.m : function f=fct(x); % fonction à intégrer f=(x-1).*exp(-x.*(x-2)); Tracé de la fonction à intégrer : >>x=0:0.01:5; >>ft=fct(x); >>plot(x,ft); >>grid; title('fonction à intégrer'); >>xlabel('x');ylabel('f(x)');  Fig. 3 : Dans le programme donné ci-dessous (quadgamm.m), on calcule l'intégrale précédente par les fonctions 'quad8' et 'gamma', pour lesquelles on compare les performances en terme de temps de calcul et du nombre d'opérations en virgule flottante. %*************************************** % Calcul de l'intégrale de f(x) par la * % la fonction 'quad8' sur un K6-233 * %*************************************** clear all; clc; flops(0); % Mise à zéro du compteur d'opérations tic; % Déclenchement du compteur temps de calcul Iquad=quad8('fct',0,5) Nbre_op1=flops, %Arrêt du compteur d’opérations tps_q=toc, %Arrêt du compteur de temps de calcul %*************************************** % Calcul de l'intégrale de f(x) par la * % la fonction 'gamma' * %*************************************** flops(0); tic; Igamma=gamma(1)/2 Nbre_op2=flops tps_g=toc En exécutant le programme quadgamm.m, on obtient les résultats suivants : >> format long >> quadgamm Iquad = 0.49999987773178 Nbre_op1 = 686 tps_q = 0.06000000000000 Igamma = 0.50000000000000 Nbre_op2 = 40 tps_g = 0.05000000000000 D'après ces résultats, on remarque bien que le nombre d'opérations en virgule flottante et le temps d'exécution dans le cas de l'utilisation de la fonction 'gamma' sont nettement inférieurs à ceux que l'on obtient en utilisant la fonction 'quad8'. La fonction 'gamma' prend des valeurs très grandes au voisinage de zéro et des entiers négatifs. Afin d'éviter un dépassement de capacité, MATLAB dispose de la fonction 'gammaln' qui retourne le logarithme népérien de ces valeurs. Exemple : >>format >>x=5e-200; >>gamma(x) ans = 2.0000e+199 >>gammaln(x) ans = 458.9076  Résolution numérique des équations différentielles et des équations aux dérivées partielles 1. Introduction Le comportement dynamique des systèmes est un sujet très important en physique. Un système mécanique par exemple, se traduit par des déplacements, des vitesses et des accélérations. Un système électrique ou électronique, se traduit par des tensions, des intensités et des dérivées temporelles sur ces quantités. Un système thermique se traduit par des températures, des gradients de température, de coefficients d’échanges thermiques, etc. En général, les équations utilisées pour décrire de tels comportements dynamiques, incluent des quantités inconnues représentant les fonctions recherchées et leurs dérivées. Une équation qui comporte une ou plusieurs dérivées de la fonction inconnue est appelée ‘équation différentielle’, qui est représentée dans MATLAB par l’abréviation ‘ODE’. L’ordre de cette équation est déterminé par l’ordre du degré le plus élevé de la dérivation. Les équations différentielles peuvent être classée en deux catégories : les équations différentielles avec des conditions initiales et les équations différentielles avec des conditions aux limites. 2. Équations différentielles du 1er ordre Les équations du 1er ordre à conditions initiales peuvent être écrites sous la forme : où est une fonction de et de , et la seconde équation représente la condition initiale sans laquelle la solution de l’équation différentielle ne peut être évaluée. Exemple 1 : Soit à résoudre l’équation différentielle suivante : où : Déterminer numériquement. On choisit un pas de temps , et on donne comme condition initiale. Tracer la solution de cette équation différentielle pour. Solution : L’équation est équivalente à où et. La solution de cette équation peut être effectuée par la méthode d’Euler, qui consiste à écrire : Le programme MATLAB suivant (eul1.m) permet la résolution de cette équation : %******************************************* % Résolution d'une équation différentielle * % par la méthode d'Euler à l'ordre 1 * %******************************************* clear ; clc ; t=0 ;n=0 ;V=0 ; C=0.27 ;M=70 ;g=9.81 ;h=0.1 ; t_R(1)=t ;V_R(1)=V ; while t