UNIVERSITE D’ABOMEY-CALAVI

--------------
ECOLE POLYTECHNIQUE D’ABOMEY-CALAVI
--------------

Introduction à la Méthode
des Eléments Finis
EDO 1D

Thomas Brice GBAGUIDI

1
Présentation générale de la Méthode des
Eléments Finis
2.1 Introduction
Pour analyser un phénomène naturel en générale ou un problème d’ingénierie
en particulier, on est souvent amené à développer un modèle mathématique
pouvant décrire d’une manière aussi fiable que possible le problème en
question.
Le développement d’un modèle mathématique s’appuie généralement sur
quelques postulats de base et plusieurs hypothèses simplificatrices pour
aboutir à des équations gouvernantes qui sont souvent des équations
différentielles auxquelles sont ajoutées des conditions aux limites. Exemple,
la théorie d’élasticité s’appuie sur le postula fondamental de l’existence du
vecteur contrainte et les équations générales d’élasticité linéaire isotrope sont
obtenues avec les hypothèses de petites déformations, d’homogénéité et
d’isotropie des matériaux ainsi que la linéarité des relations liants les
contraintes et les déformations.
La résolution analytique d’équations différentielles pose parfois des difficultés
insurmontables, et une solution exacte décrivant bien le problème étudié n’est
pas toujours facile à trouver. Le recours aux modèles physiques et à la
simulation expérimentale pour la recherche d’une solution analogue à la
solution recherchée peut s’avérer coûteux en temps et en moyens.
Avec les progrès enregistrés dans le domaine de l’informatique et les
performances des ordinateurs de plus en plus grandes, il est devenu possible
de résoudre des systèmes d’équations différentielles très complexes. Plusieurs
techniques de résolution numérique ont été ainsi développées et appliquées
avec succès pour avoir des solutions satisfaisantes à des problèmes
d’ingénierie très variés.
La méthode des éléments finis est l’une des techniques numériques les plus
puissantes. Elle est basée sur une idée simple : subdiviser (discrétiser) une
forme complexe en un grand nombre de sous-domaines élémentaires de forme
géométrique simple (éléments finis) interconnectés en des points appelés
nœuds. On retrouve les premières applications véritables de la méthode des

2
éléments finis en 1956 en mécanique des structures. Un groupe de chercheurs
(Turner, Clough, Martin et Topp [59]) de Boeing utilisent cette méthode pour
calculer la voilure d’un avion.
La méthode des éléments finis est maintenant reconnue comme l’une des
principales méthodes de résolution des équations aux dérivées partielles (EDP)
dans des géométries quelconques, que ce soit en dimension un, deux ou trois.
L’un des avantages majeurs de cette méthode est le fait qu’elle offre la
possibilité de développer un programme permettant de résoudre, avec peu de
modifications, plusieurs types de problèmes. En particulier, toute forme
complexe d’un domaine géométrique où un problème est bien posé avec toutes
les conditions aux limites, peut être facilement traité par la méthode des
éléments finis.
Cette méthode consiste à diviser le domaine physique à traiter en plusieurs
sous domaines appelés éléments finis à dimensions non infinitésimales. La
solution recherchée est remplacée dans chaque élément par une
approximation avec des polynômes simples et le domaine peut ensuite être
reconstitué avec l’assemblage ou sommation de tous les éléments.

Problème physique

Discrétisation du domaine Modèle mathématique

Choix des fonctions Formulation

d’interpolation variationnelle

Ecriture des matrices

élémentaires

Assemblage des matrices

globales, application des
conditions aux limites et
résolution du système

Figure 1.1 : étapes de la mise en œuvre de la méthode des éléments finis

3
Etape 1 : Formulation des équations gouvernantes et des conditions aux
limites.
La majorité des problèmes d'ingénierie sont décrits par des équations
différentielles aux dérivées partielles associées à des conditions aux limites
définies sur un domaine et son contour. L'application de la MEF exige une
réécriture de ces équations sous forme intégrale. La formulation faible est
souvent utilisée pour inclure les conditions aux limites.
Etape 2 : Division du domaine en sous domaines.
Cette étape consiste à discrétiser le domaine en éléments et calculer les
connectivités de chacun ainsi que les coordonnées de ses nœuds. Elle
constitue ainsi la phase de préparation des données géométriques.
Etape 3 : Approximation sur un élément.
Dans chaque élément la variable tel que le déplacement, la pression, la
température, est approximée par une simple fonction linéaire, polynomiale ou
autres. Le degré du polynôme d'interpolation est relié au nombre de nœuds de
l'élément. L'approximation nodale est appropriée. C'est dans cette étape que
se fait la construction des matrices élémentaires.
Etape 4 : Assemblage et application des conditions aux limites.
Toutes les propriétés de l'élément (masse, rigidité,...) doivent être assemblées
afin de former le système algébrique pour les valeurs nodales des variables
physiques. C'est à ce niveau qu'on utilise les connectivités calculées à l'étape
2 pour construire les matrices globales à partir des matrices élémentaires.
Etape 5 : Résolution du système global :
Le système global peut être linéaire ou non linéaire. Il définit soit un problème
d'équilibre qui concerne un cas stationnaire ou statique ou un problème de
valeurs critiques où il faut déterminer les valeurs et vecteurs propres du système
qui correspondent généralement aux fréquences et modes propres d'un système
physique.
Un problème de propagation qui concerne le cas transitoire (non stationnaire) dans
lequel il faut déterminer les variations dans le temps des variables physiques et la
propagation d'une valeur initiale. Les méthodes d'intégration pas à pas sont les plus
fréquentes telles que, méthode des différences finies centrées, méthode de Newmark,
méthode de Wilson.

4
A ces méthodes doivent être associées des techniques d'itération pour traiter le cas
non linéaire. La plus célèbre est la méthode de Newton Raphson.

2.4 Formulation variationnelle

Actuellement, le principe des travaux virtuels est bien connu et très répandu. Il est
souvent formulé en termes d'égalité des travaux effectués par les forces extérieures
et intérieures lors d'un déplacement virtuel quelconque. Ce concept est essentiel pour
la résolution des équations aux dérivées partielles. En effet, les déplacements sont
remplacés par une fonction arbitraire continue sur le domaine et l'équation est
réécrite sous forme intégrale
2.4.1 Forme forte
Un problème classique d’équations différentielles gouvernant un système
physique s'énonce comme suit:
Trouver une fonction u  V; V espace des fonctions, telle que :
A(u) = 0 |  ; B(u) = 0 |  (2.1)
où A(u) est l’ensemble d'équations gouvernantes définies sur le domaine  et
B(u) est l'ensemble des conditions aux limites que les fonctions u doivent
vérifier sur le contour  (fig.2.2). La fonction u peut être un scalaire tel que la
température, la pression, … ou un vecteur tel que le déplacement, la vitesse,
…

𝛀 B(u)=0

A(u)=0

Figure 1.2 : domaine géométrique et frontière

Le problème variationnel associé au système (2.1) s'écrit en prenant l’intégrale du
système d’équations gouvernantes pondérées par des fonctions poids, l’énoncé
devient :
Trouver u V tel que :

∀𝑤 ∈ 𝑉, ∫ w. A(𝑢)dΩ = 0 (2.2)
Ω

Cette équation est appelée forme intégrale forte de l'équation différentielle (ou du
système d’équations différentielles). Elle est analogue à l'expression des travaux
virtuels. En fait la solution de (2.2) a encore plus de portée, on peut affirmer que si

5
elle est satisfaite pour toute fonction poids w, alors l'équation différentielle (2.1) est
satisfaite en tout point du domaine .

(i) Exemple
On considère l’équation différentielle du second ordre suivante :
𝑑2 𝑢
𝐴(𝑢) = + 1 − 𝑥 = 0; 0 ≤ 𝑥 ≤ 1 (2.3)
𝑑𝑥 2
Définie dans le domaine unidimensionnel Ω = [0, 𝐿] avec les conditions aux
limites :
𝑢(𝑥 = 0) = 0 𝑒𝑡 𝑢( 𝑥 = 1 ) = 0 (2.4)
Dans ce cas 𝐵(𝑢) est l’ensemble des valeurs imposées aux deux bords du
domaine. En unidimensionnel,  se réduit à deux points. La forme intégrale
associée à l’équation A(u) s’écrit:
𝑑2 𝑢 1
𝑑2 𝑢 1
∫ 𝑤( + 1 − 𝑥) 𝑑Ω = 0; ∫ 𝑤 𝑑𝑥 = ∫ 𝑤(𝑥 − 1)𝑑𝑥 (2.5)
Ω 𝑑𝑥2 0 𝑑𝑥2 0

Avec la forme des termes à intégrer, on voit que la recherche de solutions

approximatives pour la fonction inconnue u, requiert l’utilisation de
polynômes ou de fonctions d’interpolation dérivables au moins deux fois. De
plus les conditions aux limites doivent être vérifiées à part puisqu’elles
n’interviennent pas dans l’équation intégrale ci-dessus, d’où l’introduction de
la forme intégrale faible.
2.4.2 Forme faible
Pour satisfaire les conditions aux limites nous avons deux manières de
procéder; soit par le choix de la fonction de pondération, soit vérifier que :

∫ 𝑤𝐵(𝑢) 𝑑Γ = 0 (2.6)
Γ

Dans la pratique, il est possible d'intégrer (2.2) par partie et de la remplacer par :

∫ 𝐶(𝑤)𝐷(𝑢) 𝑑Ω + ∫ 𝐸(𝑤)𝐹(𝑢) 𝑑Γ = 0 (2.7)

Ω Γ

Les opérateurs C, D, E et F contiennent des dérivées d'ordre moins élevé, d'où

un choix de fonctions d’approximation de u plus large.
Cette équation est la formulation faible de l'équation différentielle, elle forme
la base de l'approximation par éléments finis.
Remarque : Pour obtenir de telles formulations intégrales, nous disposons de
deux procédés: le premier est la méthode des résidus pondérés connue sous

6
le nom de la méthode de Galerkin, le deuxième est la détermination de
fonctionnelles variationnelles que l'on cherche à rendre stationnaire.
Dans la pratique, il n'est pas toujours facile de trouver une fonctionnelle. Le
premier procédé est plus utilisé ; il consiste à choisir 𝑤 = 𝛿𝑢 fonction de Dirac
(une perturbation de la fonction cherchée) et d'utiliser l'approximation nodale
pour la discrétisation. w s'appelle aussi fonction poids d'où le mot : "pondéré".
(i) Exemple
Pour obtenir la forme variationnelle faible de l’exemple précédent (équation
2.5) on intègre par partie le premier terme.
1 1 1
𝑑𝑤 𝑑𝑢 𝑑𝑢
−∫ 𝑑𝑥 + [𝑤 ] = ∫ 𝑤(𝑥 − 1)𝑑𝑥 (2.8)
0 𝑑𝑥 𝑑𝑥 𝑑𝑥 0 0

On voit maintenant que les conditions aux limites notamment sur les dérivées
sont explicitement prises en compte dans l’équation. Les valeurs imposées à
la fonction elle-même seront traitées lors de la résolution des systèmes
discrets.

2.5 Discrétisation du domaine

La méthode des éléments finis est une méthode d'approximation par sous
domaines, donc avant toute application il faut diviser le domaine à étudier en
éléments. Chaque élément est défini géométriquement par un nombre de
nœuds bien déterminé qui constituent en général ses sommets. (Fig. 1.3)

Frontière
Nœud interne

Elément

Nœud
frontière

Figure 1.3 : discrétisation du domaine-éléments triangulaires

La discrétisation géométrique doit respecter les règles suivantes :
• Un nœud d'un élément ne doit pas être intérieur à un côté d'un autre
du même type. (Fig. 1.4 a)

7
 • Aucun élément bidimensionnel ne doit être plat, éviter les angles
proches de 180° ou de 0°. (Fig. 1.4 b)
• Deux éléments distincts ne peuvent avoir en commun que des points
situés dans leurs frontières communes ; le recouvrement est exclu. (Fig.
1.4 c)
• L'ensemble de tous éléments doit constituer un domaine aussi proche
que possible du domaine donné ; les trous entre éléments sont exclus.
(Fig. 1.4 d)

Figure 1.4 : règle de discrétisation

Le résultat du procédé de discrétisation doit contenir deux données
essentielles qui sont les coordonnées des nœuds et les connectivités des
éléments. On doit numéroter tous les nœuds et les éléments de façon à avoir
des matrices globales à petite largeur de bande, pour cela, la numérotation se
fait selon la plus petite largeur du domaine.
2.6 Approximation sur l'élément
Après avoir défini l'élément, on peut remplacer la fonction exacte par une
approximative.
On utilise souvent des polynômes ou des fonctions faciles à mettre en œuvre
sur ordinateur.
2.6.1 Approximation polynomiale et approximation nodale
La fonction approchée est exprimée, dans le cas unidimensionnel, par :
𝑛

𝑢 = ∑ 𝑎𝑖 𝑥 𝑖 (2.9)
𝑖=0

Qu’on peut écrire sous la forme matricielle suivante :

8
 𝑎0
𝑎
𝑢 =< 1 𝑥 𝑥2 ⋯ > [𝑎1 ] ≡< 𝑃(𝑥) > {𝐴} (2.10)
2
⋮

Cette forme d'approximation est appelée interpolation polynomiale. Si on exprime la

fonction sur tous les nœuds on obtient pour chaque nœud i de coordonnée xi :

𝑢𝑖 =< 𝑃(𝑥𝑖 ) > {𝐴} ≡ ∑ 𝑃𝑖𝑗 𝑎𝑗 (2.11)

𝑗

Soit pour tous les nœuds :

< 𝑃1𝑗 > {𝑎𝑗 }
{𝑈} = [< ⋯ > {⋯ }] ≡ 𝑈𝑛 = 𝑃𝑛 𝑎𝑛 (2.12)
< 𝑃𝑛𝑗 > {𝑎𝑗 }
Un : représente les valeurs aux nœuds de la fonction.
Pn : valeurs des polynômes aux nœuds de coordonnées xi.
an : variables généralisées qui sont les facteurs des polynômes.
L’inconvénient de cette interpolation réside dans l’utilisation des paramètres
ai comme variable de base, des calculs supplémentaires sont nécessaires pour
calculer la fonction recherchée u. Afin d’éviter ces calculs, on peut mettre les
valeurs de la fonction u aux nœuds comme variables de base en procédant
comme suit :
A partir de l’équation (2.12), on peut tirer les an en fonction des un et on les
remplace dans l’équation (2.10). Ce qui donne :
𝑢 =< 𝑃(𝑥) > 𝑃𝑛−1 𝑈𝑛 ≡< 𝑁(𝑥) > 𝑈𝑛 (2.13)
C'est la forme la plus utilisée par le fait que ses variables sont les valeurs de
la fonction aux nœuds, la résolution donne directement ces valeurs.
Ce type d’approximation est appelée interpolation nodale, les fonctions Ni sont
appelées fonction de forme, elles sont fonction du type d’élément utilisé pour
la discrétisation géométrique.

9
3. Equation différentielles du premier ordre
3.2 Etapes à suivre
Dans le cas des équations différentielles, le problème est formulé directement
sous forme mathématique et revient à déterminer une fonction inconnue u
définie dans un domaine  et régie par une équation différentielle avec des
conditions aux limites
3.2.1 Formulation du problème
On prend comme exemple l’équation différentielle ordinaire d’ordre un
suivante :
𝑑𝑢
+ 2𝑥(𝑢 − 1) = 0 𝑎𝑣𝑒𝑐 𝑢(0) = 0
𝑑𝑥
Dans ce problème,  =   est un domaine de dimension 1 ; sa frontière  se
réduit à deux points : 0 et 2.
2
La solution exacte de cette équation est : 𝑢𝑒 = 1 − 𝑒 −𝑥
3.2.2 Discrétisation du domaine
Le domaine  est divisé en n segments (appelés éléments) de taille 1/n.
Chaque élément contient deux nœuds sur lesquelles la fonction u est
interpolée.

La division du domaine  en plusieurs éléments est appelée maillage. On

utilise deux tableaux pour la description du maillage : tableau de connectivités
des éléments et tableau des coordonnées des nœuds. Pour un exemple de trois
éléments on obtient les deux tableaux comme suit :
Tableau des connectivités Tableau des coordonnées
Elément Nœud 1 (début) Nœud 2 (fin) Nœud Coordonnée (x)
1 1 2 1 0.0
2 2 3 2 0.5
3 3 4 3 1.0
4 4 5 4 1.5
5 2.0

Les lignes de commandes MATLAB qui permettent d’obtenir les deux tableaux
sont :

10
n = 4; % nombre d’éléments
x = 0:2/n:2; % coordonnées des noeuds
for i = 1:n % début de boucle sur les éléments
t(i,:) = [i, 1+i]; % connectivite de chaque élément
end; % fin de boucle

ou bien, sous forme plus compacte :

x = [0:2/n:2]'
t([1:n],:) = [[1:n]',[2:n+1]']

Remarque
Pour un nombre plus élevé d’éléments il suffit d’augmenter la valeur de n .
x(i) donne la coordonnée du nœud i ; x(2) affiche : 0.5
t(i,:) donne les deux nœuds dé l’élément i ; t(2,:) affiche : 2 3
t(:,i) donne la colonne i du tableau t ; t(:,2) affiche :
2
3
4
5
3.2.3 Discrétisation et interpolation sur l’élément
On peut interpoler la fonction u recherchée dans un élément par un polynôme.
L’ordre du polynôme conditionne la précision de la solution approchée. Pour
un élément à deux nœuds on peut prendre :
𝑢 = 𝑎0 + 𝑎1 𝑥 (3.1)

soit sous forme vectorielle :

𝑎0
𝑢 =< 1 𝑥 > [𝑎 ] ≡ 𝑢 = 𝑝𝑎𝑛
1

avec p vecteur ligne contenant les monômes xn

et an vecteur colonne contenant les facteurs du polynôme.

11
Cette approximation de la fonction inconnue u est appelée interpolation polynomiale,
elle est fonction de  et  qui sont des coefficients sans valeur physique. Pour utiliser
les valeurs de u aux nœuds on cherche une interpolation en fonction de u1 et u2.
L’interpolation polynomiale aux nœuds s’écrit :
𝑢1 1 𝑥1 𝑎0
[𝑢 ] = [ ] [ ] ≡ 𝑢𝑛 = 𝑃𝑛 𝑎𝑛
2 1 𝑥2 𝑎1
L’inverse de ce système d’équations donne les paramètres an.
𝑎0 1 𝑥 −𝑥1 𝑢1
𝑎𝑛 = 𝑃𝑛−1 𝑢𝑛 ≡ [𝑎 ] = [ 2 ][ ]
1 𝑥2 − 𝑥1 −1 1 𝑢2
En remplaçant les an on peut maintenant approcher la fonction u par :
1 𝑥 −𝑥1 𝑢1 𝑥2 − 𝑥 𝑥 − 𝑥1 𝑢
𝑢 =< 1 𝑥> [ 2 ] [𝑢 ] =< > [𝑢1 ]
𝑥2 − 𝑥1 −1 1 2 𝑥2 − 𝑥1 𝑥2 − 𝑥1 2

𝑢 = 𝑁𝑈𝑛
avec N est un vecteur ligne contenant des fonctions de x appelées fonctions de
forme.
Cette interpolation est appelée interpolation nodale puisqu’elle dépend des
valeurs aux nœuds de la fonction inconnue u.
3.2.4 Propriétés des fonctions de forme
Il est intéressant de relever les propriétés suivantes pour les fonctions de forme
N:

Remarque
Les deux fonctions de forme peuvent s’écrire sous forme des polynômes de
Jacobi :
𝑥 − 𝑥𝑗
𝑁𝑖 (𝑥) = ∏ (3.2)
𝑥𝑖 − 𝑥𝑗
𝑗≠𝑖

3.2.5. Elément de type Lagrange

Ainsi les éléments pour lesquels les fonctions de forme sont des polynômes de
Jacobi sont appelés éléments de type Jacobi. On peut construire directement,

12
sans passer par l’interpolation polynomiale, des fonctions de forme pour des
éléments d’ordre supérieur (plus de deux nœuds).
Par exemple les fonctions de forme de l’élément à trois nœuds sont :
(𝑥 − 𝑥2 )(𝑥 − 𝑥3 )
𝑁1 (𝑥) =
(𝑥1 − 𝑥2 )(𝑥1 − 𝑥3 )
(𝑥 − 𝑥1 )(𝑥 − 𝑥3 )
𝑁2 (𝑥) = (3.3)
(𝑥2 − 𝑥1 )(𝑥2 − 𝑥3 )
(𝑥 − 𝑥1 )(𝑥 − 𝑥2 )
𝑁3 (𝑥) =
{ (𝑥3 − 𝑥1 )(𝑥3 − 𝑥2 )
On peut vérifier aisément sur la figure ci-dessous les propriétés aux nœuds des trois
fonctions ainsi que leur somme sur l’élément.

1 𝑠𝑖 𝑥 = 𝑥1 0 𝑠𝑖 𝑥 = 𝑥1 0 𝑠𝑖 𝑥 = 𝑥1
𝑁1 (𝑥) ={ 0 𝑠𝑖 𝑥 = 𝑥2 ; 𝑁2 (𝑥) ={ 1 𝑠𝑖 𝑥 = 𝑥2 ; 𝑁3 (𝑥) = { 𝑠𝑖 𝑥 = 𝑥2
0
0 𝑠𝑖 𝑥 = 𝑥3 0 𝑠𝑖 𝑥 = 𝑥3 1 𝑠𝑖 𝑥 = 𝑥3
Le scripte MATLAB qui permet de calculer et de tracer ces trois fonctions est le
suivant :
x1 = 2;
x2 = 5;
x3 = 8;
x = x1:0.1:x3;
N1 = (x-x2).*(x-x3)/(x1-x2)/(x1-x3);
N2 = (x-x1).*(x-x3)/(x2-x1)/(x2-x3);
N3 = (x-x1).*(x-x2)/(x3-x1)/(x3-x2);
S = N1+N2+N3;
plot(x,N1,x,N2,x,N3,x,S)
Noter qu’il est préférable d’utiliser un fichier fonction (M-file) d’usage plus
général.
3.2.5.1 Fonction m-file définissant Ni pour des points xi donnés
function [ y ] = FoncForm( i,X,x )
%=================================================
% Forme de forme pour élément de type Lagrange
%-------------------------------------------------
% Utilisation:
% [ y ] = FoncForm( i,X,x )
%-------------------------------------------------
% Description:
% Arguments d'entrée:
% i : Numéro d'ordre de la fonction de Jacobi
% X : Vecteur (ligne ou colonne) définissant

13
% les points xi
% x : Scalaire où on souhaite évaluer Ni
%
% Argument de sortie:
% y : image de x par Ni
%--------------------------------------------------
%==================================================
n=max(size(X,1),size(X,2));
p=1;
for j=1:n
if j~=i
p=p*(x-X(j))/(X(i)-X(j));
end
end
y=p;
end

Création de l’ensemble des fonctions de forme vérifiant et de leurs dérivées par un même
m-file
On sait que
𝑁𝑋
𝑥 − 𝑋𝑗
𝑁𝑖 (𝑋, 𝑥) = ∏
𝑋𝑖 − 𝑋𝑗
𝑗≠𝑖
𝑗=1

𝑁𝑋 𝑁𝑋 𝑁𝑋
𝑥 − 𝑋𝑗 𝑥 − 𝑋𝑗 𝑥 − 𝑋𝑗
⋯ ∏ ⋯
𝑁(𝑋, 𝑥) = [∏ 𝑋 − 𝑋 𝑋𝑖 − 𝑋𝑗 ∏
𝑋𝑁𝑋 − 𝑋𝑗 ]
𝑗≠11 𝑗 𝑗≠𝑖 𝑗≠𝑁𝑋
𝑗=1
𝑗=1 𝑗=1

Donc

𝑁𝑋 𝑁𝑋
𝜕 1 𝜕
[𝑁 (𝑋, 𝑥)] = ∏ ∏(𝑥 − 𝑋𝑗 )
𝜕𝑥 𝑖 𝑋𝑖 − 𝑋𝑗 𝜕𝑥
𝑗≠𝑖 𝑗≠𝑖
𝑗=1 𝑗=1
( ) ( )

𝑁𝑋 𝑁𝑋 𝑁𝑋
1
= ∏ ∑ ∏(𝑥 − 𝑋𝑗 )
𝑋𝑖 − 𝑋𝑗
𝑗≠𝑖 𝑗≠𝑖 𝑘≠𝑖
𝑗=1 𝑗=1 𝑘≠𝑗
( ) ( 𝑘=1 )

𝑁𝑋 𝑁𝑋 𝑁𝑋
1
= ∑ (∏ ) ∏(𝑥 − 𝑋𝑗 )
𝑋𝑖 − 𝑋𝑘
𝑗≠𝑖 𝑘≠𝑖 𝑘≠𝑖
𝑗=1 𝑘=1 𝑘≠𝑗
( 𝑘=1 )

14
 𝑁𝑋 𝑁𝑋 𝑁𝑋
1 1
=∑ ∏ ∏(𝑥 − 𝑋𝑗 )
𝑋𝑖 − 𝑋𝑗 𝑋𝑖 − 𝑋𝑘
𝑗≠𝑖 𝑘≠𝑖 𝑘≠𝑖
𝑗=1 𝑘≠𝑗 𝑘≠𝑗
( ( 𝑘=1 ) ( 𝑘=1 ))

𝑁𝑋 𝑁𝑋
1 𝑥 − 𝑋𝑗
=∑ ∏
𝑋𝑖 − 𝑋𝑗 𝑋𝑖 − 𝑋𝑘
𝑗≠𝑖 𝑘≠𝑖
𝑗=1 𝑘≠𝑗
( ( 𝑘=1 ))
𝑁𝑋
𝜕 1
[𝑁𝑖 (𝑋, 𝑥)] = ∑ ( 𝑁 (𝑋𝑗 , 𝑥)) (3.3)
𝜕𝑥 𝑋𝑖 − 𝑋𝑗 𝑖
𝑗≠𝑖
𝑗=1

Où 𝑋𝑗 est le vecteur obtenu en supprimant la composante j de 𝑋

𝜕
[𝑁(𝑋, 𝑥)]
𝜕𝑥
𝑁𝑋 𝑁𝑋 𝑁𝑋
1 1 1
𝑗 ⋯ ∑( 𝑁 (𝑋𝑗 , 𝑥)) ⋯ 𝑁 (𝑋𝑗 , 𝑥))]
= [∑ (𝑋 − 𝑋 𝑁1 (𝑋 , 𝑥)) 𝑋𝑖 − 𝑋𝑗 𝑖 ∑ (
𝑋𝑁𝑋 − 𝑋𝑗 𝑁𝑋
𝑗≠1 1 𝑗 𝑗≠𝑖 𝑗≠𝑁𝑋
𝑗=1
𝑗=1 𝑗=1

Création de 𝑁(𝑋, 𝑥) dans MATLAB

function [ y ] = LagMEF( X,t )
%=========================================================
% Fonctions de forme de type Jacobi
% Auteur: GBAGUIDI Thomas Brice
%=========================================================
N=length(X);
y1=zeros(N,1);
for i=1:N
y1(i,1)=1;
for k=1:N
if k~=i
y1(i,1)=y1(i,1)*(t-X(k))/(X(i)-X(k));
end
end
end
y=y1';
end

𝑑𝑁(𝑋,𝑥)
Création de 𝑑𝑥
function [y] = dLagMEF(X,t)
%=========================================================
% Dérivée première des fonctions de forme de type Jacobi
% Auteur: GBAGUIDI Thomas Brice
%=========================================================
N=length(X);

15
y1=zeros(N,1);
for i=1:N
y1(i,1)=0;
for j=1:N
if j<i
Xj=X;
Xj(j)=[];
xx=LagMEF(Xj,t);
y1(i,1)=y1(i,1)+xx(i-1)/(X(i)-X(j));
elseif j>i
Xj=X;
Xj(j)=[];
xx=LagMEF(Xj,t);
y1(i,1)=y1(i,1)+xx(i)/(X(i)-X(j));
end
end
end
y=y1';
end

Script représentant dans un même graphique les fonctions 𝑁𝑖 (𝑥) et leur somme, pour un vecteur
𝑋 = [𝑥1 , 𝑥2 , … , 𝑥𝑛 ] donné.
X=[0 0.1 0.2 0.3 0.4];
nn=size(X,2);
xx=linspace(X(1),X(end),1001);
kk=size(xx,2);
yy=zeros(nn,kk);
for k=1:kk
for ii=1:nn
yy(ii,k)=FoncForm(ii,X,xx(k));
end
end
plot(xx,yy,xx,sum(yy(1:nn,:)));
grid on;
legend('N1','N2','N3','N4','N5','Somme Ni');

On obtient le résultat suivant

1.2
N1
N2
1
N3
N4
0.8
N5
Somme Ni
0.6

0.4

0.2

-0.2

-0.4

-0.6
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

16
3.2.5.2 Construction des fonctions de forme polynomiales avec la dérivée imposée à un
nœud 𝑋𝑘 donné (interpolation de Hermite) :
𝑇
On impose par exemple que pour un vecteur de nœuds 𝑋 = [𝑋1 , 𝑋2 , … , 𝑋𝑁𝑋 ] ,
𝑁𝑖 (𝑋, 𝑘, 𝑥𝑗 ) = 𝛿𝑖𝑗
𝜕𝑁𝑖 (𝑋, 𝑘, 𝑥)
| = 0, 𝑖≠𝑘
𝜕𝑥 𝑥=𝑋 (3.4)
𝑘
𝜕𝑁𝑘 (𝑋, 𝑘, 𝑥)
| = 𝑑𝑢𝑘 ,
{ 𝜕𝑥 𝑥=𝑋
𝑘

Où 𝑘, et 𝑑𝑢𝑘 sont donnés

On a
𝑁𝑖 (𝑋, 𝑘, 𝑥𝑗 ) = 𝛿𝑖𝑗 ⇒ 𝑁𝑖 (𝑋, 𝑘, 𝑥) = 𝐴𝑖𝑘 (𝑥)𝐿𝑖 (𝑋, 𝑥) (3.5)
Cherchons les polynômes 𝐴𝑖𝑘 (𝑥) de degré 1.
𝑁𝑖 (𝑋, 𝑘, 𝑥𝑗 ) = (𝑎𝑖𝑘 𝑥𝑗 + 𝑏𝑖𝑘 )𝛿𝑖𝑗 = 𝛿𝑖𝑗 (3.6)
Alors
𝑎𝑖𝑘 𝑥𝑖 + 𝑏𝑖𝑘 = 1 (3.7)
𝜕𝑁𝑖 (𝑋, 𝑘, 𝑥) 𝜕𝐿𝑘 (𝑋, 𝑥)
= 𝑎𝑖𝑘 𝐿𝑖 (𝑋, 𝑥) + (𝑎𝑖𝑘 𝑥 + 𝑏𝑖𝑘 )
𝜕𝑥 𝜕𝑥
Ainsi
𝜕𝑁𝑖 (𝑋, 𝑘, 𝑥) 𝜕𝐿𝑖 (𝑋, 𝑥)
| = 𝛿𝑖𝑘 𝑑𝑢𝑘 ⇒ 𝑎𝑖𝑘 𝐿𝑖 (𝑋, 𝑋𝑘 ) + (𝑎𝑖𝑘 𝑋𝑘 + 𝑏𝑖𝑘 ) | = 𝛿𝑖𝑘 𝑑𝑢𝑘
𝜕𝑥 𝑥=𝑋
𝜕𝑥 𝑥=𝑋
𝑘 𝑘

D’après l’équation (3.3), on a

𝑁𝑋
𝜕𝐿𝑖 (𝑋, 𝑥) 1
= ∑( 𝐿 (𝑋𝑗 , 𝑥))
𝜕𝑥 𝑋𝑖 − 𝑋𝑗 𝑖
𝑗≠𝑖
𝑗=1

donc
𝑁𝑋
1
𝑎𝑖𝑘 𝐿𝑖 (𝑋, 𝑋𝑘 ) + (𝑎𝑖𝑘 𝑋𝑘 + 𝑏𝑖𝑘 ) ∑ ( 𝐿 (𝑋𝑗 , 𝑋𝑘 )) = 𝛿𝑖𝑘 𝑑𝑢𝑘
𝑋𝑖 − 𝑋𝑗 𝑖
𝑗≠𝑖
𝑗=1

Ou encore
𝑁𝑋 𝑁𝑋
1 1
𝑎𝑖𝑘 𝛿𝑖𝑘 + 𝑋𝑘 ∑ ( 𝐿𝑖 (𝑋𝑗 , 𝑋𝑘 )) + 𝑏𝑖𝑘 ∑ ( 𝐿 (𝑋𝑗 , 𝑋𝑘 )) = 𝛿𝑖𝑘 𝑑𝑢𝑘 (3.8)
𝑋𝑖 − 𝑋𝑗 𝑋𝑖 − 𝑋𝑗 𝑖
𝑗≠𝑖 𝑗≠𝑖
[ 𝑗=1 ] 𝑗=1

Et par substitution de (3.7) dans (3.8), on obtient

17
 𝑁𝑋 𝑁𝑋
1 1
𝑎𝑖𝑘 𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 ) ∑ ( 𝐿𝑖 (𝑋𝑗 , 𝑋𝑘 )) = 𝛿𝑖𝑘 𝑑𝑢𝑘 − ∑ ( 𝐿 (𝑋𝑗 , 𝑋𝑘 )) (3.9)
𝑋𝑖 − 𝑋𝑗 𝑋𝑖 − 𝑋𝑗 𝑖
𝑗≠𝑖 𝑗≠𝑖
[ 𝑗=1 ] 𝑗=1

Pour 𝒊 ≠ 𝒌, on a
L’égalité (3.8) devient après simplification des termes nuls,
1 1
𝑎𝑖𝑘 [(𝑋𝑘 − 𝑋𝑖 ) 𝐿𝑖 (𝑋 𝑘 , 𝑋𝑘 )] = − 𝐿 (𝑋 𝑘 , 𝑋𝑘 ) (3.10)
𝑋𝑖 − 𝑋𝑘 𝑋𝑖 − 𝑋𝑘 𝑖
On sait que 𝑥 ⟼ 𝐿𝑖 (𝑋 𝑘 , 𝑥) a exactement n -1 racines réelles distinctes qui sont les composantes
de 𝑋 𝑘 , n étant le nombre de composante de 𝑋 𝑘 . Mais, 𝑋𝑘 n’est pas une racine de 𝐿𝑖 (𝑋 𝑘 , 𝑥) car,
c’est la composante de X supprimée pour constituer 𝑋 𝑘 (voir paragraphe 3.2.5.1). On peut donc
simplifier l’égalité (3.10), ce qui donne :
1
𝑎𝑖𝑘 =
𝑋𝑖 − 𝑋𝑘
−𝑋𝑘
𝑏𝑖𝑘 = 1 − 𝑎𝑖𝑘 𝑋𝑖 =
𝑋𝑖 − 𝑋𝑘
𝑥 − 𝑋𝑘
𝑁𝑖 (𝑋, 𝑘, 𝑥) = 𝐿 (𝑋, 𝑥) (3.11)
𝑋𝑖 − 𝑋𝑘 𝑖
Pour 𝒊 = 𝒌, :
L’égalité (3.9) devient
𝑁𝑋 𝑁𝑋
1 1
𝑎𝑘𝑘 1 + (𝑋𝑘 − 𝑋𝑘 ) ∑ ( 𝐿𝑘 (𝑋𝑗 , 𝑋𝑘 )) = 𝑑𝑢𝑘 − ∑ ( 𝐿 (𝑋𝑗 , 𝑋𝑘 )) (3.12)
𝑋𝑘 − 𝑋𝑗 𝑋𝑘 − 𝑋𝑗 𝑘
𝑗≠𝑘 𝑗≠𝑘
[ 𝑗=1 ] 𝑗=1

𝐿𝑘 (𝑋𝑗 , 𝑋𝑘 ) = 1 ∀𝑗 ≠ 𝑘
Donc
𝑁𝑋
1
𝑎𝑘𝑘 = 𝑑𝑢𝑘 − ∑ ( )
𝑋𝑘 − 𝑋𝑗
𝑗≠𝑘
𝑗=1

𝑁𝑋
1
𝑏𝑘𝑘 = 1 − 𝑋𝑘 𝑑𝑢𝑘 − ∑ ( )
𝑋𝑘 − 𝑋𝑗
𝑗≠𝑖
[ 𝑗=1 ]

𝑁𝑋
1
𝑁𝑘 (𝑋, 𝑘, 𝑥) = 𝑑𝑢𝑘 − ∑ ( ) (𝑥 − 𝑋𝑘 ) + 1 𝐿𝑘 (𝑋, 𝑥) (3.13)
𝑋𝑘 − 𝑋𝑗
𝑗≠𝑖
[( 𝑗=1
) ]

18
M-file
function [ y ] = Ndu0( X,du0k,k,t )
%=========================================================
% Fonctions de forme de type polinomiale Ni
% tel que
% ° Ni(X,k,xj)=𝛿ij
% ° dNi(X,k,xj)/dx=0 si i≠k
% ° dNk(X,k,xj)/dx=duok
%
% Auteur: GBAGUIDI Thomas Brice
%=========================================================
N=length(X);
y1=zeros(N,1);
% y2=y1;
y2=0;
for j=1:N
if j~=k
y2=y2+1/(X(k)-X(j));
end
end
for i=1:N
if i==k
y1(i,1)=(du0k-y2)*(t-X(k))+1;
else
y1(i,1)=(t-X(k))/(X(i)-X(k));
end
for j=1:N
if j~=i
y1(i,1)=y1(i,1)*(t-X(j))/(X(i)-X(j));
end
end
end
y=y1';
end

Exemple de programme pour la représentation graphique des Ni

clear all;
X=[0 0.1 0.2 0.3];
t=linspace(X(0),X(end),1001);
n=[size(X,2),size(t,2)];
y=zeros(n(1),n(2));
for i=1:n(2)
y(:,i)=Ndu0( X,-10,4,t(i) );
end
plot(t,y);
grid on;

On obtient successivement pour 𝑘 = 4 avec 𝑑𝑢0𝑘 = 5, puis 𝑘 = 2 avec 𝑑𝑢0𝑘 = 2

19
1.5
N1 est de dérivée 1.5
première égale à N2 est de dérivée
5 en x0=0 première égale
1 1 à 2 en x1 = 0.25

0.5 0.5

0 0

-0.5 -0.5
N1 N1, N3, N4 et N5 sont
N2, N3, N4 et N5 N1
N2 de dérivées premières
sont de dérivées nulles en x1=0.25 N2
-1 N3 -1
premières nulles N3
en x0=0 N4
N4
N5
N5
-1.5 -1.5
0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Les dérivées
Pour 𝒊 ≠ 𝒌
𝑑 1 𝑥 − 𝑋𝑘 𝑑
𝑁𝑖 (𝑋, 𝑘, 𝑥) = 𝐿𝑖 (𝑋, 𝑥) + 𝐿 (𝑋, 𝑥)
𝑑𝑥 𝑋𝑖 − 𝑋𝑘 𝑋𝑖 − 𝑋𝑘 𝑑𝑥 𝑖
𝑁𝑋
𝑑 1 𝑥 − 𝑋𝑘 1
𝑁𝑖 (𝑋, 𝑘, 𝑥) = 𝐿𝑖 (𝑋, 𝑥) + ∑( 𝐿𝑖 (𝑋𝑗 , 𝑥))
𝑑𝑥 𝑋𝑖 − 𝑋𝑘 𝑋𝑖 − 𝑋𝑘 𝑋𝑖 − 𝑋𝑗
𝑗≠𝑖
𝑗=1

Pour 𝒊 = 𝒌

𝑁𝑋
𝑑 𝑑 𝑑 1
𝑁𝑖 (𝑋, 𝑘, 𝑥) = 𝑁𝑘 (𝑋, 𝑘, 𝑥) = 𝑑𝑢𝑘 − ∑ ( ) (𝑥 − 𝑋𝑘 ) + 1 𝐿𝑘 (𝑋, 𝑥)
𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑋𝑘 − 𝑋𝑗
𝑗≠𝑖
[( 𝑗=1
) ]

𝑁𝑋
𝑑 1
𝑁 (𝑋, 𝑘, 𝑥) = 𝑑𝑢𝑘 − ∑ ( ) 𝐿𝑘 (𝑋, 𝑥)
𝑑𝑥 𝑘 𝑋𝑘 − 𝑋𝑗
𝑗≠𝑖
𝑗=1
( )

𝑁𝑋 𝑁𝑋
1 1
+ 𝑑𝑢𝑘 − ∑ ( ) (𝑥 − 𝑋𝑘 ) + 1 ∑ ( 𝐿 (𝑋𝑗 , 𝑥)) (3.14)
𝑋𝑘 − 𝑋𝑗 𝑋𝑘 − 𝑋𝑗 𝑘
𝑗≠𝑖 𝑗≠𝑘
[( 𝑗=1
) ] 𝑗=1

20
 25

20 dN1/dx
dN2/dx
15 dN3/dx
dN4/dx
10 dN5/dx

-5

-10

-15

-20

-25
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

1.5

0.5

-0.5

N1
-1 N2
N3
N4
N5
-1.5
0 0.2 0.4 0.6 0.8 1

Généralisation
On pose ici
𝑁𝑖 (𝑋, 𝑘, 𝑥𝑗 ) = 𝛿𝑖𝑗
{ 𝜕𝑁𝑖 (𝑋, 𝑘, 𝑥) (3.15)
| = 𝑑𝑢𝑖 ,
𝜕𝑥 𝑥=𝑋 𝑘

Avec
𝑁𝑖 (𝑋, 𝑘, 𝑥) = (𝑎𝑖𝑘 𝑥 + 𝑏𝑖𝑘 )𝐿𝑖 (𝑋, 𝑥)
𝜕𝐿𝑖 (𝑋,𝑥)
En posant | = 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ), on obtient le système suivant
𝜕𝑥 𝑥=𝑋𝑘

𝑎𝑖𝑘 [𝛿𝑖𝑘 + 𝑋𝑘 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )] + 𝑏𝑖𝑘 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) = 𝑑𝑢𝑖

{
𝑎𝑖𝑘 𝑋𝑖 + 𝑏𝑖𝑘 = 1
Ce qui donne

𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )

𝑎𝑖𝑘 =
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )

21
 𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝛿𝑖𝑘 + 𝑋𝑘 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) − 𝑋𝑖 𝑑𝑢𝑖
𝑏𝑖𝑘 = 1 − 𝑋𝑖 =
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )

𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝛿𝑖𝑘 + 𝑋𝑘 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) − 𝑋𝑖 𝑑𝑢𝑖

𝑁𝑖 (𝑋, 𝑘, 𝑥) = [ 𝑥+ ] 𝐿 (𝑋, 𝑥)
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑖

𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )

𝑁𝑖 (𝑋, 𝑘, 𝑥) = [ (𝑥 − 𝑋𝑘 )
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝑢𝑖
+ ] 𝐿 (𝑋, 𝑥) (3.16)
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑖

𝜕𝑁𝑖 (𝑋, 𝑘, 𝑥) 𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )

= 𝐿 (𝑋, 𝑥)
𝜕𝑥 𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑖
𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
+[ (𝑥 − 𝑋𝑘 )
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝑢𝑖
+ ] 𝑑𝐿𝑖 (𝑋, 𝑥) (3.17)
𝛿𝑖𝑘 + (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )

Pour 𝒊 ≠ 𝒌, 𝜹𝒊𝒌 = 𝟎 et on a
𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
𝑎𝑖𝑘 =
(𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
𝑋𝑘 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) − 𝑋𝑖 𝑑𝑢𝑖
𝑏𝑖𝑘 =
(𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
On obtient
𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑑𝑢𝑖
𝑁𝑖 (𝑋, 𝑘, 𝑥) = [ (𝑥 − 𝑋𝑘 ) + ] 𝐿 (𝑋, 𝑥) (3.18)
(𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑖
Avec le cas particulier 𝑑𝑢𝑖 = 0 , on retrouve la formule (3.11).
𝜕𝑁𝑖 (𝑋, 𝑘, 𝑥) 𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
= 𝐿 (𝑋, 𝑥)
𝜕𝑥 (𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑖
𝑑𝑢𝑖 − 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑑𝑢𝑖
+[ (𝑥 − 𝑋𝑘 ) + ] 𝑑𝐿𝑖 (𝑋, 𝑥) (3.19)
(𝑋𝑘 − 𝑋𝑖 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 ) 𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )

Pour 𝒊 = 𝒌, 𝜹𝒊𝒌 = 𝟏 et on a
𝑑𝑢𝑘 − 𝑑𝐿𝑘 (𝑋, 𝑋𝑘 )
𝑎𝑘𝑘 = = 𝑑𝑢𝑘 − 𝑑𝐿𝑘 (𝑋, 𝑋𝑘 )
1 + (𝑋𝑘 − 𝑋𝑘 )𝑑𝐿𝑘 (𝑋, 𝑋𝑘 )
1 + 𝑋𝑘 𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ) − 𝑋𝑘 𝑑𝑢𝑘
𝑏𝑘𝑘 = = 1 + 𝑋𝑘 (𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ) − 𝑑𝑢𝑘 )
1 + (𝑋𝑘 − 𝑋𝑘 )𝑑𝐿𝑖 (𝑋, 𝑋𝑘 )
On obtient

22
 𝑁𝑘 (𝑋, 𝑘, 𝑥) = [(𝑑𝑢𝑘 − 𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ))𝑥 + 1 + 𝑋𝑘 (𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ) − 𝑑𝑢𝑘 )]𝐿𝑘 (𝑋, 𝑥)
𝑁𝑘 (𝑋, 𝑘, 𝑥) = [(𝑑𝑢𝑘 − 𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ))(𝑥 − 𝑋𝑘 ) + 1]𝐿𝑘 (𝑋, 𝑥) (3.20)

𝜕𝑁𝑘 (𝑋, 𝑘, 𝑥)
= (𝑑𝑢𝑘 − 𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ))𝐿𝑘 (𝑋, 𝑥)
𝜕𝑥
+[(𝑑𝑢𝑘 − 𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ))(𝑥 − 𝑋𝑘 ) + 1]𝑑𝐿𝑘 (𝑋, 𝑥) (3.21)

Les formules (3.19) et (3.20) sont respectivement équivalente à (3.13) et (3.14).

𝑁𝑋
1
𝑑𝐿𝑘 (𝑋, 𝑋𝑘 ) = ∑
𝑋𝑘 − 𝑋𝑗
𝑗=1
𝑗≠𝑘

3.2.6 Matrices élémentaires (interpolation linéaire)

Le calcul des matrices élémentaires passe par la réécriture du problème sous forme intégrale :
𝑑𝑢
∫ 𝛿𝑢 [ + 2𝑥(𝑢 − 1)] 𝑑Ω
Ω 𝑑𝑥
avec 𝛿𝑢 une fonction de pondération prise égale à une perturbation de la fonction inconnue u.
Le domaine  comprend l’intervalle 0 à 2, 𝑑Ω = 𝑑𝑥 et avec l’interpolation nodale on a :
𝑑𝑢 𝑑𝑁
= 𝑈 ; 𝑒𝑡 𝛿𝑢 = 𝑁 𝛿𝑈𝑛
𝑑𝑥 𝑑𝑥 𝑛
Puisque seules les fonctions N dépendent de x et les perturbations ne touchent que les valeurs
de u, pour commodité on écrit : 𝛿𝑢 = (𝛿𝑈𝑛 )𝑇 𝑁 𝑇
L’intégrale de 0 à 2 peut être remplacée par la somme des intégrales de xi à xi+1 (ou bien :
l’intégrale sur  est la somme des intégrales sur e, avec e le domaine de chaque élément)
2 𝑥𝑖+1 𝑥2
∫ = ∑ ∫ ≡∫ = ∑ ∫ = ∑ ∫
Ω é𝑙é𝑚𝑒𝑛𝑡 Ωe 0 𝑖=1,𝑛 x𝑖 é𝑙é𝑚𝑒𝑛𝑡 x1

La forme intégrale de l’équation différentielle devient alors pour chaque élément :

𝑥2 𝑥2 𝑥2
𝑑𝑢
∫ [𝛿𝑢 ] 𝑑𝑥 + ∫ [𝛿𝑢 2𝑥 𝑢]𝑑𝑥 − ∫ [𝛿𝑢] 2𝑥 𝑑𝑥 = 0
𝑑𝑥
𝑥1 𝑥1 𝑥1
𝑥2 𝑥2 𝑥2
𝑑𝑁
∫ [𝛿𝑈𝑛𝑇 𝑁 𝑇 𝑈 ] 𝑑𝑥 + ∫ [𝛿𝑈𝑛𝑇 𝑁 𝑇 2𝑥 𝑁 𝑈𝑛 ]𝑑𝑥 − ∫ [𝛿𝑈𝑛𝑇 𝑁 𝑇 ] 2𝑥 𝑑𝑥 = 0
𝑑𝑥 𝑛
𝑥1 𝑥1 𝑥1

Cette écriture discrétisée est valable pour tous les types d’éléments. Dans le cas particulier d’un
élément linéaire à deux nœuds, elle s’écrit comme suit :
𝑥2
𝑁1 𝑢
< 𝛿𝑢1 𝛿𝑢2 > { ∫ [ ] < 𝑑𝑁1 𝑑𝑁2 > 𝑑𝑥} [𝑢1 ] +
𝑁2 2
𝑥1

23
 𝑥2
𝑁1 𝑢
< 𝛿𝑢1 𝛿𝑢2 > { ∫ [ ] 2𝑥 < 𝑁1 𝑁2 > 𝑑𝑥} [𝑢1 ] −
𝑁2 2
𝑥1
𝑥2
𝑁1
< 𝛿𝑢1 𝛿𝑢2 > { ∫ [ ] 2𝑥𝑑𝑥} = 0
𝑁2
𝑥1

On voit qu’il est possible de simplifier (𝛿𝑢𝑛 )𝑇 puisqu’il n’est pas nul et revient à chaque terme.
Au fait c’est les autres termes qui doivent être nuls ! on le voit bien avec cet exemple :
𝑣1 𝑤1
< 𝛿𝑢1 𝛿𝑢2 > [𝑣 ] +< 𝛿𝑢1 𝛿𝑢2 > [𝑤 ] = 0
2 2

donne: 𝛿𝑢1 𝑣1 + 𝛿𝑢2 𝑣2 + 𝛿𝑢1 𝑤1 + 𝛿𝑢2 𝑤2 ; soit : 𝑣1 + 𝑤1 = 0 et 𝑣2 + 𝑤2 = 0

𝑣1 𝑤1
ou : [𝑣 ] + [𝑤 ] = 0
2 2

Finalement l’équation intégrale discrétisée se met sous la forme matricielle :

𝑥2 𝑥2 𝑥2
𝑁 𝑁 𝑢 𝑁
{ ∫ [ 1 ] < 𝑑𝑁1 𝑑𝑁2 > 𝑑𝑥 + ∫ [ 1 ] 2𝑥 < 𝑁1 𝑁2 > 𝑑𝑥} [𝑢1 ] − ∫ [ 1 ] 2𝑥𝑑𝑥 = 0
𝑁2 𝑁2 2 𝑁2
𝑥1 𝑥1 𝑥1
𝑥2 𝑥2 𝑥2
𝑁 𝑑𝑁 𝑁1 𝑑𝑁2 𝑁 2𝑥𝑁1 𝑁1 2𝑥𝑁2 𝑢1 2𝑥𝑁1
{∫ [ 1 1 ] 𝑑𝑥 + ∫ [ 1 ] 𝑑𝑥 } [𝑢 ] = ∫ [ ] 𝑑𝑥
𝑁2 𝑑𝑁1 𝑁2 𝑑𝑁2 𝑁2 2𝑥𝑁1 𝑁2 2𝑥𝑁2 2 2𝑥𝑁2
𝑥1 𝑥1 𝑥1

Qui est un système d’équations linéaires 𝐾𝑒𝑈𝑒 = 𝐹𝑒 ; avec Ke et Fe sont appelés matrice et
vecteur élémentaires du système d’équation. Dans le cas de la présente équation différentielle
K est la somme de deux matrices : Ke = Ke1 + Ke2 tel que :
𝑥2 𝑥2
𝑁 𝑑𝑁 𝑁1 𝑑𝑁2 𝑁 2𝑥𝑁1 𝑁1 2𝑥𝑁2
𝐾𝑒1 = ∫ [ 1 1 ] 𝑑𝑥 ; 𝐾𝑒2 = ∫ [ 1 ] 𝑑𝑥
𝑁2 𝑑𝑁1 𝑁2 𝑑𝑁2 𝑁2 2𝑥𝑁1 𝑁2 2𝑥𝑁2
𝑥1 𝑥1

En remplaçant les fonctions de forme et leurs dérivées par leurs expressions respectives on
obtient :
𝑥2𝑥 − 𝑥2 1 𝑥 − 𝑥2 1
𝑥 − 𝑥2 𝑥1 − 𝑥2 𝑥1 − 𝑥2 𝑥2 − 𝑥1
𝐾𝑒1 = ∫ 1 𝑑𝑥
𝑥 − 𝑥1 1 𝑥 − 𝑥1 1
𝑥1
[𝑥2 − 𝑥1 𝑥1 − 𝑥2 𝑥2 − 𝑥1 𝑥2 − 𝑥1 ]
𝑥2 𝑥 − 𝑥2 𝑥 − 𝑥2 𝑥 − 𝑥2 𝑥 − 𝑥1
2𝑥 2𝑥
𝑥 −𝑥 𝑥1 − 𝑥2 𝑥1 − 𝑥2 𝑥2 − 𝑥1
𝐾𝑒2 = ∫ [ 𝑥1 − 𝑥 2 𝑥 − 𝑥 𝑥 − 𝑥 𝑥 − 𝑥1 ] 𝑑𝑥
1 2 1
𝑥1 2𝑥 2𝑥
𝑥2 − 𝑥1 𝑥1 − 𝑥2 𝑥2 − 𝑥1 𝑥2 − 𝑥1
Soit après intégration des composantes des deux matrices :
1 −1 1
𝐾𝑒1 = [ ]
2 −1 1

24
 𝑥2 − 𝑥1 3𝑥1 + 𝑥2 𝑥1 + 𝑥2
𝐾𝑒2 = [ ]
6 𝑥1 + 𝑥2 𝑥1 + 3𝑥2
Le vecteur F est donné par :
𝑥2
2𝑥(𝑥 − 𝑥2 )/(𝑥1 − 𝑥2 )
𝐹𝑒 = ∫ [ ] 𝑑𝑥
2𝑥(𝑥 − 𝑥1 )/(𝑥2 − 𝑥1 )
𝑥1

soit :
𝑥2 − 𝑥1 2𝑥1 + 𝑥2
𝐹𝑒 = [ ]
3 𝑥1 + 2𝑥2
3.2.7 Remarque
Il est possible d’utiliser MATLAB pour intégrer analytiquement les matrices élémentaires, Le script
peut être le suivant :
syms x x1 x2 real % déclaration de variables symboliques
N = [(x-x2)/(x1-x2) (x-x1)/(x2-x1)] % fonctions de forme
dN = simple(diff(N,x)) % dérivées des fonctions de forme
Ke1 = simplify( int(N' * dN , x, x1, x2) ) % matrice Ke1
Ke2 = simplify( int(N' * 2*x * N , x, x1, x2) ) % matrice Ke2
Fe = simplify( int(N' * 2*x , x, x1, x2) ) % vecteur Fe

3.3 Etape 4 : Assemblage

Le calcul des matrices élémentaires permet d’obtenir pour les quatre éléments les systèmes
d’équations élémentaires suivants :
Elément 1 : x = 0 ; x = ½ ;

(1) − 1⁄2 1⁄ 1 1
2] ; 𝐾𝑒 (1) = ⁄2 − 0 [3 × 0 + ⁄2 0 + 1⁄2
𝐾𝑒1 = [ 2 ]
− 1⁄2 1⁄
2
6 0 + 1⁄2 0 + 3 × 1⁄2

(1) (1) 1 −11 13 −0.4583 0.5417

𝐾𝑒 (1) = 𝐾𝑒1 + 𝐾𝑒2 = [ ]=[ ]
24 −11 15 −0.4583 0.6250
1⁄ − 0 2 × 0 + 1⁄
𝐹𝑒 (1) = 2 [ 2] = 1 [1] = [0.0833]
3 0 + 2 × 1⁄2 12 2 0.1667

1 −11 13 𝑢1 1 1
𝐾𝑒 (1) 𝑈𝑒 (1) = 𝐹𝑒 (1) ; [ ] [𝑢 ] = [ ]
24 −11 15 2 12 2
On notant par U les valeurs de la fonction u aux cinq nœuds, les valeurs du vecteur élémentaire
𝑈𝑒 (1) de l’élément 1 correspondent aux composantes 𝑢1 et 𝑢2 vecteur global U, celles de
𝑈𝑒 (2) de l’élément 2 correspondent à la seconde et troisième du vecteur global U, celles de
𝑈𝑒 (3) à la troisième et quatrième et celles de 𝑈𝑒 (4) à la quatrième et cinquième.
Elément 2 : x1 = 1/2 ; x2 = 1 ;
1 −7 15 −0.2917 0.6250
𝐾𝑒 (2) = [ ]=[ ]
24 −9 19 −0.3750 0.7917

25
 1 4 0.3333
𝐹𝑒 (2) = [ ]=[ ]
12 5 0.4167
1 −7 15 𝑢1 1 4
𝐾𝑒 (2) 𝑈𝑒 (2) = 𝐹𝑒 (2) ; [ ] [𝑢 ] = [ ]
24 −9 19 2 12 5
Elément 3 : x1 = ; x2 = 3/2 ;
1 −3 17 −0.1250 0.7083
𝐾𝑒 (3) = [ ]=[ ]
24 −7 23 −0.2917 0.9583
1 7 0.5833
𝐹𝑒 (3) = [ ]=[ ]
12 8 0.6667
1 −3 17 𝑢1 1 7
𝐾𝑒 (3) 𝑈𝑒 (3) = 𝐹𝑒 (3) ; [ ] [𝑢 ] = [ ]
24 −7 23 2 12 8
Elément 4 : x1= 3/2 ; x2 = ;
1 1 19 −0.0417 0.7917
𝐾𝑒 (4) = [ ]=[ ]
24 −5 27 −0.2083 1.1250
1 10 0.8333
𝐹𝑒 (4) = [ ]=[ ]
12 11 0.9167
1 1 19 𝑢1 1 10
𝐾𝑒 (4) 𝑈𝑒 (4) = 𝐹𝑒 (4) ; [ ] [𝑢 ] = [ ]
24 −5 27 2 12 11
En réécrivant les systèmes élémentaires en fonction de toutes les composantes de U on obtient
:
Elément 1 :
−11 13 0 0 0 𝑢1 1
1 −11 15 0 0 0 𝑢 2 1 2
0 0 0 0 0 𝑢3 = 0
24 0 0 0 0 0 4 𝑢 12 0
[ 0 0 0 0 0] [𝑢5 ] [0]
Elément 2 :
0 0 0 0 0 𝑢1 0
1 0 −7 15 0 0 𝑢2 1 4
0 −9 19 0 0 𝑢3 = 5
24 0 0 0 0 0 4 𝑢 12 0
[0 0 0 0 0] [𝑢5 ] [0]
Elément 3 :
0 0 0 0 0 𝑢1 0
1 0 0 0 0 0 𝑢2 1 0
0 0 −3 17 0 𝑢3 = 7
24 0 0 −7 23 0 4 𝑢 12 8
[0 0 0 0 0] [𝑢5 ] [0]
Elément 4 :

26
 0 0 0 0 0 𝑢1 0
1 0 0 0 0 0 𝑢2 1 0
0 0 0 0 0 𝑢3 = 0
24 0 0 0 1 19 4 𝑢 12 10
[0 0 0 −5 27] [𝑢5 ] [11]
En prenant maintenant la somme ( somme des intégrales), le système global s’écrit enfin :

−11 13 0 0 0 𝑢1 1
1 −11 15 − 7 15 0 0 𝑢2 1 2+4
0 −9 19 − 3 17 0 𝑢3 = 5 + 7 ≡ 𝐾𝑈 = 𝐹
24 0 0 −7 23 + 1 19 𝑢4 12 8 + 10
[ 0 0 0 −5 27] [𝑢5 ] [ 11 ]
On appel cette opération assemblage. En pratique, on ne procède pas de cette manière pour des
raisons d’économie de mémoire et de temps de calcul mais on fait un assemblage des matrices
élémentaires en utilisant les connectivités des éléments. La matrice globale K est d’abord
initialisée à la matrice nulle, ensuite à chaque construction de matrice élémentaire, on localise
avec une table là où il faut l’ajouter à la matrice globale. Dans le cas d’un degré de liberté par
nœud (ce cas présent) la table de localisation correspond à la table des connectivités. Exemple
; pour l’élément 1 la table de localisation est t = 1 2 on ajoute donc la matrice 𝐾𝑒 (1) à K(1,1)
, K(1,2) , K(2,1) et K(2,2) ; et pour l’élément 2 la table de localisation est t = 2 3 on ajoute la
matrice 𝐾𝑒 (2) à K(2,2) , K(2,3) , K(3,2) et K(3,3). Et ainsi de suite.
Le script MATLAB qui permet de faire ce type d’assemblage est simple :

K = zeros(n+1); % initialisation de la matrice globale (n+1 noeuds)

F = zeros(n+1,1) ; % initialisation du vecteur global (remarquer ,1)
for i = 1:n % boucle sur les éléments
t = [i i+1]; % table de localisation de chaque élément
x1 = x(i); x2 = x(i+1); % coordonnées des deux noeuds de l’élément
[Ke,Fe] = MatElt2Nd(x1,x2); %Fonction pour calculer la mat. et vect. élementaires
K(t,t) = K(t,t) + Ke; % Assemblage de la matrice globale
F(t) = F(t) + Fe; % Assemblage du vecteur global
end; % fin de boucle

Programme utilisant integrale, LagMEF( X,t ) et dLagMEF( X,t )

%==========================================================================
% Programme pour la résolution des équations de
% la forme du/dx+a(x)u+b(x)=0 avec u(x0)=u0
% par la méthode des éléments finis sur
% l'intervalle [x0,xn]
% Fonctions de forme: polynôme d'interpolation
% de Lagrange
% Auteur: GBAGUIDI Thomas Brice
%==========================================================================
clc;
clear all;
% Intervalle
x0=-2;xn=2;
% Nombre de subdivision de l'intervalle [x0,xn] pour la MEF

27
n=20;
u0=2;
% Nombre de subdivision de l'intervalle [x0,xn] pour les intégrales
Nint=1000;
% Subdivision
X=linspace(x0,xn,n+1);
% Degré polynôme
nn=4; % Nombre de noeuds par élément
% Les fonctions a(t) et b(t)
a=@(t)(2*t); % Définition de la fonction a(t)
b=@(t)(0); % Définition de la fonction b(t)
% La matrice Ke et le vecteur Fe
Ke=zeros(n+1,n+1); % Initialisation de Ke à 0
Fe=zeros(n+1,1); % Initialisation de Fe à 0
XX=zeros(nn+1,n-nn+1);
for i=1:n-nn+1
for jj=1:nn+1
XX(jj,i)=X(i+jj-1);
end
ke=@(t)(LagMEF(XX(:,i),t)'*(dLagMEF(XX(:,i),t)+a(t)*LagMEF(XX(:,i),t)));
fe=@(t)(LagMEF(XX(:,i),t)'*b(t));
intke=integrale(ke,XX(1,i),XX(nn+1,i),Nint);
intfe=integrale(fe,XX(1,i),XX(nn+1,i),Nint);
Ke(i:i+nn,i:i+nn)=Ke(i:i+nn,i:i+nn)+intke; % Assemblage de Ke
Fe(i:i+nn,1)=Fe(i:i+nn,1)+intfe; % Assemblage de Fe
end
K=Ke(2:n+1,2:n+1);%exclusion de la première colonne et première ligne de Ke
F=Fe(2:n+1)-u0*Ke(2:n+1,1);% prise en compte de la condition aux limites
u=K\F;
X1=linspace(x0,xn,1001);y3=2*exp(-X1.^2+4);
f=@(s1,s2)[-a(s1)*s2+b(s1)];
[x2,y2]=RK4expSys(f,x0,xn,u0,15);
u=[u0;u];% Prise en compte de u(x0)=u0
plot(X1,y3,'r',X,u,'ro',x2,y2,'b+');%Représentation graphique
text(0,0,'du/dx+2x*u=0,x?[0,2],u(0)=2');
legend('Solution analytique','Solution MEF','Solution RK4');
grid on ;

function y = integrale( f,a,b,n )

%=======================================================================
% Fonction pour le calcul d'intégrale par la méthode des trapèzes
% Entrées :
% f: fonction à intégrer
% a, b: borne inférieure et borne supérieure respectivement
% n: nombre de subdivision de l'intervalle d'intégration
% Sortie : l'intégrale calculée
% Auteur : GBAGUIDI Brice
%=======================================================================
h=(b-a)/n;
x0=f(a);
s=x0;
x0=f(b);
s=(s+x0)/2;
x0=a;
for i=1:n-1
x0=x0+h;
s=s+f(x0);
end
y=s*h;
end

28
 120

du/dx+2x*u=0, Solution analytique

sur [0,2],avec u(0)=2
Solution MEF
Solution RK4

100

80

60

40

20

0
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

3.4 Etape 5a : Résolution – Application des CAL

Avant de résoudre le système il faut appliquer les conditions aux limites de la fonction u.
Au nœud 1 de coordonnée x = 0, u = 0 ; ce qui se traduit par u = 0. Ceci induit la réduction du nombre
d’équations total à résoudre. Le système global devient alors :
15 − 7 15 0 0 𝑢2 2+4
1 −9 19 − 3 17 0 𝑢3 1 5+7
[ ][ ] = [ ]
24 0 −7 23 + 1 19 𝑢4 12 8 + 10
0 0 −5 27 𝑢5 11
La ligne et la colonne d’indice 1 qui correspondent à la valeur u1 ont été supprimées de la matrice K et
du vecteur F. Si 𝑢(0) = 𝑎 ≠ 0 (𝑢1 = 𝑎) alors on retranche de F le produit de la 1ère colonne de K par
a.
La commande MATLAB qui permet de supprimer une ligne ou une colonne d’une matrice consiste à
lui affecter un vecteur nul (vide) :
A(i, :) = [] % supprime la ligne i de la matrice A
A(:, i) = [] % supprime la colonne i
V(i) =[] % supprime la composante i du vecteur V

Pour appliquer la condition aux limites u = 0, il suffit d’écrire :

K(1,:) = []; % suppression de la 1ere ligne de K

K(:,1) = []; % suppression de la 1ere colonne de K
F(1,:) = []; % suppression de la 1ere ligne du vecteur F

Remarque : pour u(0) = a, on insère, avant les suppressions, la commande :

F = F - K(:,1)*a.

3.5 Etape 5b : Résolution – Calcul de la solution

29
La solution peut être maintenant calculée par un des algorithmes de résolution de système linéaires.
Notons que dans ce cas la matrice K est tridiagonale, mais avec MATLAB il suffit d’utiliser la division
gauche U = K \ F. Les résultats sont donnés dans le tableau et la figure suivants :

Coord. x Solution Exacte Solution MEF Erreur (%)

0.5000 0.2212 0.2488 12.4555
1.0000 0.6321 0.6673 5.5705
1.5000 0.8946 0.9154 2.3225
2.0000 0.9817 0.9843 0.2694

0.8
Solution exacte
0.6 Solution MEF

0.4

0.2

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

3.6 Etude de la convergence

Toute étude qui passe par une solution numérique approchée doit faire l’objet d’un examen de
convergence. Dans les études par éléments finis la convergence peut être atteinte en augmentant le degré
des polynômes d’interpolation (nombre de nœuds par élément), en utilisant un nombre d’éléments assez
grand avec un raffinement du maillage ou bien en cherchant une bonne adaptation du maillage au
problème traité. La figure si après donne le pourcentage d’erreur relative au point x = 0.5 en fonction du
nombre n d’éléments utilisés. On voit que le décroissement de l’erreur relative est une fonction de type
1/n et tends à s’annuler d’une manière monotone. Un nombre d’éléments n  30 suffit pour avoir une
erreur acceptable.

function [U, Ue, Err, x] = ExempleEquaDiff(n)

%
% du/dx - 2*x*(1-u) = 0
% u(0) = 0
%
%----------------------------------------------------
% Solution avec des fonctions d'interpolation à 2Nds
%----------------------------------------------------
x = [0:2/n:2]'; % vecteur des coordonnées
K = zeros(n+1 ) ; % initialisations
F = zeros(n+1,1) ;

30
for i = 1:n % boucle sur les éléments
j = i+1;
t = [i j];
x1 = x(i);
x2 = x(j);
[Ke, Fe] = MatElt2Nd(x1,x2);
K(t,t) = K(t,t) + Ke;
F(t) = F(t) + Fe;
end;
K(1,:) = []; % application des CAL
K(:,1) = [];
F(1,:) = [];
U = K\F; % résolution
U =[0;U]; % incorporation de la valeur initiale pour plot
Ue = 1-exp(-x.^2); % calcul de la solution exacte
Err = 100*(U-Ue)./Ue; % calcul de l’erreur relative
plot(x,U,x,Ue) % trace les deux solutions
end

function [Ke, Fe] = MatElt2Nd(x1,x2)

%-------------------------------------------
% Calcul de la matrice Ke et du vecteur Fe
%-------------------------------------------
Ke1 = 1/2*[ -1 1 % matrice Ke1
-1 1 ] ;
Ke2 = (x2-x1)/6* [ 3*x1+x2 x1+x2 % matrice Ke2
x1+x2 x1+3*x2];
Ke = Ke1 + Ke2; % matrice Ke
Fe = (x2-x1)/3 * [2*x1+x2 ; x1+2*x2]; % vecteur Fe
end

La commande : [U, Ue, Err, x] = ExempleEquaDiff(30) permet d’avoir à la fois la

solution par éléments finis U, la solution exacte Ue, l’erreur relative Err et les coordonnées x avec
30 éléments.
Remarque :
Noter que ce programme permet pratiquement de résoudre toutes les équations différentielles d’ordre 1
avec la condition aux limites u(0) = 0. Seule la fonction MatElt2Nd(x1,x2) nécessite un
changement pour les matrices élémentaires.

3.7 Devoir N°2

1) Modifier le programme ExempleEquaDiff(n) pour résoudre la même équation en utilisant des
fonctions de forme quadratiques (élément à trois nœuds)
2) Ecrire une fonction MATLAB qui permet d’avoir les expressions des fonctions de forme 1D.
3) Ecrire un programme MATLAB d’éléments finis qui permet de résoudre dans le domaine [0, 3]
𝑑𝑢 2
l’équation : 𝑑𝑥 + (2𝑥 − 1)𝑢 = 0 avec (0) = 1 . Comparer avec la solution exacte 𝑢 = 𝑒 𝑥−𝑥 et
examiner la convergence au point x = 0.5.
3.8 Equations différentielles du 2nd ordre
Comme exemple d’équation du 2nd ordre, on va traiter dans le domaine  = [0 , 1], l’équation :
𝑑2 𝑢 𝑑𝑢
2
+6 + 9𝑢 = 𝑥(1 − 𝑥)
𝑑𝑥 𝑑𝑥

31
 𝑑𝑢
Avec les conditions aux limites 𝑢(0) = 0 et 𝑑𝑥 | =0
𝑥=1
1 1 −3𝑥
La solution exacte est : 𝑢 = (3𝑥 + 4)(1 − 𝑥) + 𝑒 (𝑥𝑒 3 + 8 − 12𝑥)
27 54

Le domaine  sera subdivisé (maillé) en n éléments linéaires à 2 nœuds. La forme intégrale forte de
cette équation s’écrit :
1
𝑑2 𝑢 1
𝑑𝑢 1 1
∫ 𝛿𝑢 𝑑𝑥 + 6 ∫ 𝛿𝑢 𝑑𝑥 + 9 ∫ 𝛿𝑢𝑢𝑑𝑥 = ∫ 𝛿𝑢𝑥(1 − 𝑥)𝑑𝑥
0 𝑑𝑥 2 0 𝑑𝑥 0 0

Lorsqu’il est question d’équations différentielles du second ordre et plus, on est confronté à deux
difficultés. L’interpolation des dérivées et la satisfaction des conditions aux limites. En effet plus l’ordre
des dérivées est grand, plus fort doit être le degré des polynômes ou des fonctions d’interpolation à
utiliser. Un élément linéaire à deux nœuds ne peut donc être utilisé pour ce type d’équation. De plus, la
𝑑𝑢
forme intégrale forte ne fait pas apparaître la condition 𝑑𝑥 | =0
𝑥=1

Ces deux difficultés ont conduit à l’utilisation de la forme intégrale faible qu’on peut obtenir, dans le
cas présent avec intégration par partie du premier terme :
1
𝑑2 𝑢 1
𝛿𝑑𝑢 𝑑𝑢 𝑑𝑢 1
∫ 𝛿𝑢 2 𝑑𝑥 = − ∫ 𝑑𝑥 + [𝛿𝑢 ]
0 𝑑𝑥 0 𝑑𝑥 𝑑𝑥 𝑑𝑥 0
On voit immédiatement l’avantage majeur qu’offre cette expression. Elle réduit l’ordre des dérivées
𝑑𝑢
(d’où le terme faible) et permet de prendre en compte la condition aux limite |
𝑑𝑥 𝑥=1
=0

De plus puisque 𝑢 (0) = 0, le terme 𝛿𝑢 s’annule aussi à 𝑥 = 0 (il ne peut y avoir de perturbations
dans des valeurs connues ou nulles).
Ainsi la forme intégrale s’écrit avec la prise en compte des conditions aux limites comme suit :
1 1 1 1
𝛿𝑑𝑢 𝑑𝑢 𝑑𝑢
−∫ 𝑑𝑥 + 6 ∫ 𝛿𝑢 𝑑𝑥 + 9 ∫ 𝛿𝑢𝑢𝑑𝑥 = ∫ 𝛿𝑢𝑥(1 − 𝑥)𝑑𝑥
0 𝑑𝑥 𝑑𝑥 0 𝑑𝑥 0 0

Sa discrétisation donne pour les intégrales sur les éléments :

1
𝑑𝑁 𝑇 𝑑𝑁 1
𝑑𝑁 1 1
[− ∫ 𝑑𝑥 + 6 ∫ 𝑁 𝑇 𝑑𝑥 + 9 ∫ 𝑁 𝑇 𝑁𝑑𝑥 ] 𝑈𝑛 = ∫ 𝑁 𝑇 𝑥(1 − 𝑥)𝑑𝑥
0 𝑑𝑥 𝑑𝑥 0 𝑑𝑥 0 0

Avec cette écriture il est possible d’utiliser un élément linéaire ; la matrice élémentaire Ke est la somme
de trois matrices : Ke = Ke + Ke2 + Ke3 , avec :

𝑥 − 𝑥2
𝑥2 1 𝑥 − 𝑥2 1
𝑥 − 𝑥2 𝑥1 − 𝑥2 𝑥1 − 𝑥2 𝑥2 − 𝑥1
𝐾𝑒1 = ∫ 1 𝑑𝑥
𝑥 − 𝑥1 1 𝑥 − 𝑥1 1
𝑥1
[𝑥2 − 𝑥1 𝑥1 − 𝑥2 𝑥2 − 𝑥1 𝑥2 − 𝑥1 ]

1
𝑥2 𝑥2
𝑥1 − 𝑥2 1 1 −1 1 −1
= −∫ < > 𝑑𝑥 = ∫ [ ] 𝑑𝑥
𝑥1
1 𝑥1 − 𝑥2 𝑥2 − 𝑥1 2
(𝑥1 − 𝑥2 ) 𝑥1 −1 1
[𝑥2 − 𝑥1 ]
32
 1 −1 1
𝐾𝑒1 = [ ]
𝑥2 − 𝑥1 1 −1
𝑥 − 𝑥2
𝑥2 𝑥2
𝑥 −𝑥 1 1 6 𝑥 − 𝑥2 𝑥2 − 𝑥
𝐾𝑒2 = 6 ∫ [ 𝑥1 − 𝑥 2 ] < > 𝑑𝑥 = 2
∫ [𝑥 − 𝑥 𝑥 − 𝑥1 ] 𝑑𝑥
𝑥1
1 𝑥1 − 𝑥2 𝑥2 − 𝑥1 (𝑥1 − 𝑥2 ) 𝑥1 1
𝑥2 − 𝑥1

−3 3
𝐾𝑒2 = [ ]
−3 3
𝑥 − 𝑥2
𝑥2 𝑥 − 𝑥2 𝑥 − 𝑥1
𝑥 −𝑥
𝐾𝑒3 = 9 ∫ [ 𝑥1 − 𝑥 2 ] < > 𝑑𝑥
𝑥1
1 𝑥1 − 𝑥2 𝑥2 − 𝑥1
𝑥2 − 𝑥1
𝑥2
9 (𝑥 − 𝑥2 )2 −(𝑥 − 𝑥2 )(𝑥 − 𝑥1 )
= ∫ [ ] 𝑑𝑥
(𝑥1 − 𝑥2 )2 𝑥1 −(𝑥 − 𝑥2 )(𝑥 − 𝑥1 ) (𝑥 − 𝑥1 )2

(𝑥2 − 𝑥1 )3 (𝑥1 − 𝑥2 )3
9 3 6
=
(𝑥1 − 𝑥2 ) (𝑥1 − 𝑥2 )3
2 (𝑥2 − 𝑥1 )3
[ 6 3 ]

3(𝑥2 − 𝑥1 ) 2 1
𝐾𝑒3 = [ ]
2 1 2
et le vecteur élémentaire Fe est :
𝑥 − 𝑥2
𝑥2
𝑥 −𝑥
𝐹𝑒 = ∫ [ 𝑥1 − 𝑥 2 ] . 𝑥. (1 − 𝑥)𝑑𝑥
1
𝑥1
𝑥2 − 𝑥1

𝑥1 − 𝑥2 3𝑥12 − 4𝑥1 + 2𝑥1 𝑥2 − 2𝑥2 + 𝑥22

𝐹𝑒 = [ 2 ]
12 𝑥1 − 2𝑥1 + 2𝑥1 𝑥2 − 4𝑥2 + 3𝑥22

3.9 Programme MATLAB

Remarquer les modifications introduites dans le programme précédent pour traiter cette équation.

function [U, Ue] = EquaDiff2(n)

%==================================
% d²u/dx² + 6 du/dx + 9 u = x(1-x)
% u(0) = 0 du(1)= 0
%==================================
x = [0:1/n:1]'; % modification d’1 borne d’intégration
K = zeros(n+1 ) ;
F = zeros(n+1,1) ;
for i = 1:n
j = i+1;
t = [i j];
x1 = x(i);

33
 x2 = x(j);
[Ke,Fe] = MatElt2Nd(x1,x2);
K(t,t) = K(t,t) + Ke;
F(t) = F(t) + Fe;
end;
K(1,:) = [];
K(:,1) = [];
F(1) = [];
U = K\F;
U = [0.0;U];
t = 0:0.01:1;
Ue = 1/27*(1-t).*(3*t-4)+1/54*exp(-3*t).*(8+t*exp(3)-12*t);
plot(x,U,'-.',t,Ue)
end

function [Ke,Fe] = MatElt2Nd(x1,x2)

%==========================================
% Calcul de la matrice Ke et du vecteur Fe
%==========================================
Ke1 = 1/(x2-x1)*[ -1 1 % les modifications ne touchent
1 -1 ] ; % essentiellement que les matrices
Ke2 = [ -3 3 % élémentaires
-3 3 ] ;
Ke3 = (x2-x1)* [3 3/2
3/2 3];
Ke = Ke1 + Ke2 + Ke3;
Fe = (x2-x1)/12* [(4-3*x1-2*x2)*x1+(2-x2)*x2;
(4-3*x2-2*x1)*x2+(2-x1)*x1 ];
end

Un appel de ce programme avec la commande EquaDiff2(10), puis EquaDiff2(100) permet d’avoir

les figures suivantes :
-3 -3
x 10 x 10
16 16
Solution MEF avec 10 éléments Solution MEF avec 100 éléments
14 Solution exacte 14 Solution exacte

12 12

10 10

8 8

6 6

4 4

2 2

0 0

-2 -2

-4 -4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

(a) Solution avec n=10 (b) Solution avec n=100

3.10 Comparaison avec l’intégration pas à pas

Comparativement aux méthodes d’intégration pas à pas qui traitent les problèmes aux valeurs initiales
(ou problèmes d’évolution), la méthode des éléments finis traite les problèmes aux valeurs aux limites.
Ceci se traduit par les conditions aux limites ou conditions initiales exigées par les deux méthodes.
A titre de comparaison nous allons traiter l’équation différentielle

34
 𝑑2 𝑢
+𝑢−𝑥 =0
𝑑𝑥 2
𝑑𝑢
Avec les conditions aux limites 𝑢(0) = 0 et 𝑑𝑥 | = 0 pour la MEF et les conditions initiales 𝑢(0) =
𝑥=4𝜋
𝑑2 𝑢
0 et | = 0 pour l’intégration Pas à Pas
𝑑𝑥 2 𝑥=0

La solution exacte étant ue = x − sin x. La commande MATLAB qui permet d’avoir cette solution est :
>> ue = simplify(dsolve('D2u + u -t','u(0)=0','Du(4*pi)=0'))

La fonction MATLAB d’intégration pas à pas des systèmes à un degré de liberté avec la méthode
d’accélération linéaire est utilisée pour résoudre numériquement l’équation. Un exposé détaillé de cette
méthode peut être trouvé dans l’ouvrage de dynamique des structures de Penzien et Clough.

function [u,v,a] = SdofTH(m,c,k,acc,dt)

%=======================================================
% [u,v,a] = SdofTH(m,c,k,acc,dt)
% u, v, a : histoire de la réponse en déplacements u,
% vitesses v et accélérations a d’un système à 1DDL
% m, c, k : masse, amortissement et rigidité
% acc : accélération à la base en vecteur
% dt : pas de temps
%=======================================================
u(1) = 0.0;
v(1) = 0.0;
n = length(acc);
a=zeros(1,n);
a(1) = - acc(1);
K = k + 6/dt^2 * m + 3/dt * c;
for i = 1:n-1
A = 6/dt * v(i) + 3.0 * a(i);
B = 3.0 * v(i) + dt/2 * a(i);
P = m * ( acc(i) - acc(i+1) + A ) + B * c;
du = P/K;
dv = 3/dt * du - 3.0 * v(i) - dt/2 * a(i);
u(i+1) = u(i) + du;
v(i+1) = v(i) + dv;
P = - c * v(i+1) - k * u(i+1) - m * acc(i+1);
a(i+1) = P/m;
end
end

Les commandes MATLAB qui permettent de tracer les deux solutions sont les suivantes :

clear
clc
n = 20 % nombre de pas d’intégration ou d’éléments finis
dt = 4*pi/n; % pas de temps d’intégration
t = 0:dt:4*pi; % discrétisation du temps
acc = -t; % force d’inertie Fi = - m *acc
up = sdofth(1,0,1,acc,dt); % appel de la fonction d’intégration pas à pas
[um,xm] = EquadiffMEF(n); % solution par éléments finis
x = 0:0.01:4*pi; % coordonnées x pour tracer la solution exacte
ue= x-sin(x); % solution exacte
plot(t,up,'.',xm,um,'o', x,ue) % trace les trois solutions

Remarque : La fonction EquadiffMEF(n) a la même structure que EquaDiff2(n). La différence réside

seulement dans le calcul des matrices et vecteur élémentaires.
35
Remarquer aussi les variables de sortie !
Le tracé des trois solutions avec un nombre de pas d’intégration égal au nombre d’éléments montre sur
les deux figures ci-dessous que l’intégration pas à pas est un peu plus proche de la solution exacte quand
le nombre d’éléments est petit et les deux solutions coïncident presque parfaitement avec la solution
exacte lorsque le nombre d’éléments est grand.
La différence dans le cas d’un petit nombre d’éléments est due au fait que l’élément utilisé est linéaire,
donc le polynôme d’interpolation est du première ordre. Par contre, l’intégration pas à pas utilise une
approximation avec un polynôme cubique.

10 éléments et 10 pas d’intégration 20 éléments et 20 pas d’intégration

3.11 Devoir N°3

1) Modifier le programme EquaDiff2(n) pour obtenir EquadiffMEF(n). Calculer avec MATLAB les
matrices et vecteur élémentaires et changer la fonction MatElt2Nd(x1,x2) pour calculer les nouvelles
matrices et le nouveau vecteur.
2) Introduire dans MATLAB la fonction SdofTH et retracer les deux figures précédentes
3) Utiliser la fonction FF1D(n) ci-dessous pour calculer les fonctions de forme d’un élément quadratique
(élément à trois nœuds) et calculer les expressions des matrices et vecteur élémentaires de l’équation
𝑑2 𝑢
précédente 𝑑𝑥 2 + 𝑢 − 𝑥 = 0 : avec les mêmes CAL.
function [N,dN] = FF1D(n)
%-------------------------------------------------------------------
% [N,dN] = FFDim1(n)
% N : fonction de forme d’élément unidimensionnel de type Jacobi
% dN : dérivés des fonctions de forme
% n : nombre de noeuds de l’élément
% la fonction retourne l’expression symbolique des fonctions de forme
%
%-------------------------------------------------------------------
x = sym('x','real'); % caractère x
for i=1:n % boucle pour ajouter les caractères 1 2 3 ...
c = char(48+i); % attention à n > 9
xn(i) = sym(['x',c],'real');
end;
for i=1:n % boucle pour les N(i)
N(i) = sym(1); % initialisation à 1
for j=1:n % boucle pour le produit
if j~=i
N(i) = N(i)*(x-xn(j))/(xn(i)-xn(j));
end

36
 end
end
N = expand(N); N = simplify(N); % réarrangement des termes de N
dN = simplify(diff(N,x)); % calcul des dérivée
end

3.12 Programme général pour les équations du 2nd ordre

La structure des programmes MATLAB réalisés jusqu’ici ne permet pas de traiter toutes les équations
différentielles du second ordre. Des modifications plus ou moins légères sont nécessaires pour chaque
cas, elles touchent essentiellement les matrices élémentaires et les conditions aux limites.
Il est possible de penser à mettre la fonction de calcul des matrices élémentaires dans un autre fichier et
de créer une autre fonction de traitement des conditions aux limites pour garder le reste du programme
dans une fonction qui ne sera pas modifiée. Mais, en plus du fait que le type d’élément reste fixé, les
expressions des matrices élémentaires elles-mêmes posent toujours problème puisqu’il faut les calculer
par intégration sur l’élément.
Or, dans certains cas justement, l’intégration peut être difficile ou quasiment impossible, l’intégration
numérique est une solution qui convient parfaitement notamment si on l’associe à un élément à bornes
fixes appelé élément parent ou élément de référence.
Par ailleurs, les coefficients de l’équation peuvent être eux aussi des fonctions de la variable x. Pour
mettre en œuvre un programme qui traite sans modifications des équations très diversifiées, on doit donc
en premier lieu caractériser les variables qui sont des données du problème.
Un problème aux valeurs limites défini dans un domaine Ω = [𝑥0 , 𝑥𝑙 ] et régie par une équation
différentielle linéaire du 2nd ordre s’écrit comme suit :
𝑑2 𝑢 𝑑𝑢
+ 𝑐(𝑥) + 𝑘(𝑥)𝑢 = 𝑓(𝑥)
𝑑𝑥 2 𝑑𝑥
𝑑𝑢
avec 𝑢(𝑥0 ) = 𝑢0 et | = 𝑑𝑢𝑙
𝑑𝑥 𝑥=𝑙

Si le terme du second ordre de l’équation est multiplié par une fonction m(x) alors il est préférable de
diviser toute l’équation par ce terme sinon la formulation variationnelle faible fera intervenir un terme
intégrale de plus contenant la dérivée de m(x) par rapport à x.
Les paramètres de l’équation sont les fonctions k, c et f ainsi que les valeurs 𝑢0 et 𝑑𝑢𝑙 , ceux du domaine
 (de la géométrie) sont 𝑥0 et 𝑥𝑙 et ceux propres à la méthode sont le type de l’élément caractérisé par
le nombre de nœud qu’il contient nne et le nombre d’éléments total net. Il est préférable de regrouper
tous ces paramètres dans une structure de données appelée Eq comme Equation.
Le programme d’éléments finis qu’on veut obtenir peut être appelé avec une préparation de données de
la manière suivante :
On voit clairement qu’avec ce type de script, il est facile de traiter des équations très diverses ; il suffit
de donner l’expression des coefficients dans les fonctions coef_k, coef_c et func_f, les bornes du
domaine ainsi que les conditions aux limites peuvent être changées à l’aide des valeurs de Eq.xi, Eq.xf,

37
Eq.u0 et Eq.Du. De plus nous avons maintenant le choix du type d’éléments et de leur nombre avec les
variables Eq.nne et Eq.net.
Remarque
La solution exacte de l’équation différentielle traitée comme exemple peut être obtenue avec ces
commandes :
>>ue = dsolve('D2u + 4* u = cos(t)','u(0)=1','Du(pi)=1');
>>ue = simplify(ue)
ue = 1/3*cos(t)+2/3*cos(2*t)+1/2*sin(2*t)
Pour arriver à une telle flexibilité, le programme doit être assez complet. Sa structure
générale est :
Début du programme
Géométrie (coordonnées et connectivités)
Initialisation des matrices
Début de boucle sur les éléments
Appel de fonction pour calculer  dNT dN dx
Appel de fonction pour calculer  NT c(x) dN dx
Appel de fonction pour calculer  NT k(x) N dx
Appel de fonction pour calculer  NT f(x) dx
Assemblage de la matrice K
Assemblage du vecteur F
Fin de boucle sur les éléments
Application des conditions aux limites
Résolution
Fin du programme

Le fichier fonction MefEquaDiff est le suivant :

function [U,p] = MefEquaDiff(Eq)
%============================================================
% [U,p] = MefEquaDiff(Eq)
% U : vecteur contenant les valeurs de la solution aux noeuds
% p : coordonnées des noeuds
% Eq: structure contenant les paramètres de l’équation
%
%============================================================
[t,p] = maillage1d(Eq.net,Eq.nne,Eq.xi,Eq.xf); % connectivités et
coordonnées
nnt = length(p); % nombre de noeuds total
K = zeros(nnt); % initialisation de la matrice K
F = zeros(nnt,1); % initialisation du vecteur F
for ie = 1:Eq.net % boucle sur les éléments
te = t(ie,:); % connectivité de l’élément actuel
X = p(te)'; % coordonnées de l’élément actuel
Me = -intdNTdN(X); % calcul de l’intégrale ∫ dNT. dNdx
Ce = intNTcxdN(Eq.c,X); % calcul de l’intégrale ∫ NT. c. dNdx
Ke = intNTkxN (Eq.k,X); % calcul de l’intégrale ∫ NT. k(x). Ndx
Fe = intNTfx (Eq.f,X); % calcul de l’intégrale ∫ NT. f(x)dx
K(te,te) = K(te,te) + Me + Ce + Ke; % assemblage de la matrice K
F(te ) = F(te ) + Fe ; % assemblage du vecteur F
end; % fin de boucle
F = F -K(:,1)*Eq.u0; % application des CAL
F(nnt) = F(nnt) - Eq.Du;
K(1,:) = [];
K(:,1) = [];
F(1) = [];
U = K\F; % résolution
U = [Eq.u0;U];

38
end

La première chose qu’on remarque dans cette version du programme est que le domaine est maillé avec
la fonction maillage1d qui permet de calculer à la fois les coordonnées des nœuds et les connectivités
des éléments. Cette façon de faire est particulièrement utile pour les cas de deux ou trois dimensions,
notamment lorsqu’on veut utiliser un fichier de maillage résultat d’un programme de maillage
automatique. Pour ce cas simple, on propose cette fonction comme suit :
function [t,p] = maillage1d(net,nne,xi,xf)
%=====================================================
% [t,p] = maillage1d(net,nne)
% t : table des connectivités des éléments linéiques
% p : table des coordonnées des noeuds
% net : nombre d'éléments total
% nne : nombre de noeuds par élément
% xi,xf : bornes du domaine d'intégration
%
%=====================================================
t(1,:) = 1:nne;
for i = 2:net
t(i,:) = t(i-1,:)+nne-1;
end;
n = (nne-1)*(net);
dx = (xf-xi)/n;
p = xi:dx:xf;
end
La deuxième remarque est plus importante, elle concerne le calcul des matrices élémentaires. La
fonction MatElt2N utilisée dans les versions précédentes est maintenant remplacée par quatre fonctions
dont trois : intdNTdN, intNTcxdN, intNTkxN retournent les matrices élémentaires après intégration et
la quatrième intNTfx fait le calcul du vecteur force élémentaire. Les arguments d’entrée sont un vecteur
X des coordonnées de l’élément, les pointeurs Eq.k, Eq.c et Eq.f qui prennent les adresses (ou handle)
des fonctions k(x), c(x) et f(x). La taille du vecteur X n’étant pas spécifiée, les fonctions prennent en
charge tout type d’élément linéique.
Il est à rappeler que c’est dans ces fonctions que s’effectue l’intégration des matrices, et c’est justement
ce qui a limité les versions précédentes, donc c’est ici qu’on s’attend à des modifications. En effet,
l’intégration exacte des composantes des matrices est remplacée par une somme pondérée des valeurs
de ces composantes en des points spécifiques selon la méthode d’intégration numérique de Gauss. Cette
méthode qui sera exposée plus loin, nous permet d’évaluer l’intégrale d’une fonction à l’aide d’une
somme, soit :
𝑏 𝑛

∫ 𝑓(𝑥)𝑑𝑥 = ∑ 𝑤𝑖 𝑓(𝑥𝑖 )
𝑎 𝑖=1

Avec f une fonction quelconque de x, n est le nombre de points xi dans lesquels on évalue la fonction et
𝑤𝑖 sont des valeurs de pondération. Les coordonnées xi ainsi que les poids 𝑤𝑖 sont calculés et la précision
de cette intégration numérique dépend du nombre de points n, le résultat est exacte si la fonction f est

39
un polynôme d’ordre 2n. On donne ci-dessous les points et poids pour l’intégration sur l’intervalle
−1,+1 :
n=3
x -0.34641016151378 0 0.34641016151378
w 0.55555555555556 0.88888888888889 0.55555555555556
n=4
x -0.861136311594053 0.861136311594053 -0.339981043584856
0.339981043584856
w 0.347854845137454 0.347854845137454 0.652145154862546
0.652145154862546
n=6
x -0.932469514203152 0.932469514203152 -0.661209386466265
0.661209386466265
w 0.171324492379170 0.171324492379170 0.360761573048139
0.360761573048139
x -0.238619186083197 0.238619186083197
w 0.467913934572691 0.467913934572691

La fonction MATLAB implémentée pour le programme est basée sur un nombre de point n = 6, ce qui
nous permet d’utiliser pratiquement sans perte de précision des éléments ayant jusqu’à 11 nœuds.
function [xsi,w] = wpgL6p
%====================================================================
% [xsi,w] = wpgL6p
% xsi, w : points et poids de Gauss entre -1 et +1 en vecteurs lignes
% nombre de poins : 6 points
%====================================================================
xsi = [ -0.932469514203152 , 0.932469514203152 , ...
-0.661209386466265 , 0.661209386466265 , ...
-0.238619186083197 , 0.238619186083197 , ...
];
w = [ 0.171324492379170 , 0.171324492379170 , ...
0.360761573048139 , 0.360761573048139 , ...
0.467913934572691 , 0.467913934572691 , ...
];
end

Avec cette méthode d’intégration, les expressions exactes des matrices élémentaires
ne sont plus nécessaires, d’autant plus qu’on ne connaît pas, à priori, les fonctions
à intégrer ; les coefficients de l’équation peuvent être quelconques.

Le dernier point qui reste consiste à faire le passage entre les bornes d’intégration de
l’élément réel et les bornes -1 et +1 de l’élément de référence. La transformation
géométrique :

𝑥 = 𝑁𝜉

avec  la coordonnée dans l’élément de référence, offre la meilleure solution du fait

qu’elle utilise les fonctions de forme N définies aussi dans l’élément de référence. On
peut vérifier aisément que si 𝜉 = ±1, 𝑥 = 𝑥1 , 𝑥2 pour un élément linéaire et x x 1, x 3
pour un élément quadratique, … etc.
Avec cette transformation géométrique, les dérivées des fonctions de forme dans les
deux éléments sont liées par la relation :

40
 𝑑𝑁 𝑑𝑁 𝑑𝑥
=
𝑑𝜉 𝑑𝑥 𝑑𝜉
𝑑𝑁 𝑑𝑁
Soit sous forme plus générale 𝑑𝜉
= 𝐽 𝑑𝑥 . (dans le cas de deux ou trois dimensions,
est une matrice jacobienne)

𝑑𝑥 𝑑 𝑑𝑁
𝐽= = (𝑁𝑥𝑛 ) = 𝑥
𝑑𝜉 𝑑𝜉 𝑑𝜉 𝑛
1 𝑑𝑁 1
Exemple d’un élément linéaire : 𝑁 = <1−𝜉 1+𝜉 > , = < −1 1 >
2 𝑑𝜉 2

1 𝑥1 1
et 𝐽 = 2 < −1 1 > [𝑥 ] = 2 (𝑥2 − 𝑥1 )
2

Les bornes des intégrales sont ainsi transformées = 1:

𝑥2 1
∫ 𝑓(𝑥)𝑑𝑥 = ∫ 𝑓(𝑁(𝜉))𝐽𝑑𝜉
𝑥1 −1

Finalement, avec ces deux modifications majeures qui concernent la transformation

géométrique et l’intégration numérique, les quatre fonctions des matrices et vecteur
élémentaires auront la structure suivante :
Début de fonction
Points et poids d’intégration de Gauss
Initialisation de la matrice (vecteur) élémentaire
Début de boucle sur les points de Gauss
Evaluation des fonctions de forme au point de Gauss
Calcul de J = dN * Xn
Calcul de 𝑑𝑁/𝑑𝑥 = 𝐽−1 𝑑𝑁
Somme des produits 𝑑𝑁 𝑇 . 𝑑𝑁, 𝑑𝑁 𝑇 . 𝑁, 𝑁 𝑇 . 𝑁 avec multiplication par
les poids et J
Fin de boucle sur les points de Gauss
Fin de la fonction

Le code source des quatre fonctions est le suivant :

function Me = intdNTdN(Xn)
%===================================
% integration de : dN' * dN
%===================================
[xsi, w] = wpgL4p;
npg = length(xsi);
nne = length(Xn);
Me = zeros(nne);
for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
J = dN * Xn;
dN = dN * inv(J);
Me = Me + det(J) * w(i) * (dN' * dN);
end
end

function Ce = intNTcxdN(cx,Xn)
%===================================
% integration de : N' * c(x) * dN
%===================================
[xsi, w] = wpgL4p;
npg = length(xsi);
nne = length(Xn);

41
Ce = zeros(nne);
for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
x = N * Xn;
J = dN * Xn;
dN = dN * inv(J);
c = feval(cx,x);
Ce = Ce + det(J) * w(i) * c * (N' * dN);
end
end

function Ke = intNTkxN(kx,Xn)
%==================================
% intégration de : N' * k(x) * N
%==================================
[xsi, w] = wpgL4p;
npg = length(xsi);
nne = length(Xn);
Ke = zeros(nne);
for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
x = N * Xn;
J = dN * Xn;
k = feval(kx,x);
Ke = Ke+det(J)*w(i)*k*(N'*N);
end
end

function Fe = intNTfx(fx,Xn)
%==================================
% integration de : N' * f(x)
%==================================
[xsi, w] = wpgL4p;
npg = length(xsi);
nne = length(Xn);
Fe = zeros(nne,1);
for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
x = N * Xn;
J = dN * Xn;
f = feval(fx,x);
Fe = Fe + det(J) * w(i) * f * N';
end
end

Les fonctions de forme ainsi que leurs dérivées pour des éléments de référence d’ordre
n sont :
function [N,dN] = FF1D(x,n)
%=======================================================
% [N,dN] = FF1D(x,n)
% N : fonctions de forme d'éléments linéiques 1D
% dN : dérivées des fonctions de forme
% x : point d'evaluation des fonctions
% n : nombre de noeuds de l'élément
%
%
%=======================================================

42
if isa(x,'sym') % détermine si x est symbolique
one = sym(1);
zero = sym(0);
else
one = 1.0;
zero = 0.0;
if (x < -1) || (x > 1)
error('coordonné x erronée');
end
end
xn =-1:2/(n-1):1;
for i=1:n
p(i) = one;
q(i) = one;
s(i) = zero;
for j=1:n
if j~=i
p(i) = p(i) * (x - xn(j));
q(i) = q(i) * (xn(i) - xn(j));
r = one;
for k =1:n
if (k~=i && k~=j)
r = r*(x - xn(k));
end
end
s(i) = s(i) + r;
end
end
end
N = p./q;
dN = s./q;
end

3.13 Devoir N°4

𝑑2 𝑢
1) Reprendre l’équation du devoir N°3, question 3 : 𝑑𝑥 2 + 𝑢 − 𝑥 = 0 avec ce programme.

Modifier les fonctions coef_k, coef_c et func_f, ainsi que les bornes du domaine et les
conditions aux limites. Faire une étude de convergence en type et en nombre
d’éléments.
2) Résoudre avec le programme MefEquaDiff l’équation :
𝑑2 𝑢 1 𝑑𝑢 1 + 𝑥 2
+ =
𝑑𝑥 2 𝑥 𝑑𝑥 𝑥
𝑑𝑢
avec les conditions aux limites : 𝑢(0) = 1 et 𝑑𝑥 | = −2
𝑥=4𝜋

3.14 Résolution des problèmes d’ordre 2 avec des conditions aux limites de type Robin
Il s’agit des problèmes de la forme
𝑑2 𝑢 𝑑𝑢
2
+ 𝑎(𝑥) + 𝑏(𝑥)𝑢 = 𝑐(𝑥), 𝑥 ∈ ]𝑥0 , 𝑥𝑛 [
𝑑𝑥 𝑑𝑥
𝛼1 𝑢(𝑥0 ) + 𝛽1 𝑢′(𝑥0 ) = 𝛾1 (𝛼1 , 𝛽1 ) ≠ (0,0)
{𝛼2 𝑢(𝑥𝑛 ) + 𝛽2 𝑢′(𝑥𝑛 ) = 𝛾2 (𝛼2 , 𝛽2 ) ≠ (0,0)
En suivant la démarche du paragraphe 3.8, on peut écrire :
(𝐾𝑒1 + 𝐾𝑒2 + 𝐾𝑒3 )𝑈𝑛 = 𝐹𝑒
Avec

43
 𝑥2
𝑑𝑁 𝑇 𝑑𝑁
𝐾𝑒1 = − ∫ 𝑑𝑥
𝑥1 𝑑𝑥 𝑑𝑥
𝑥2
𝑑𝑁
𝐾𝑒2 = ∫ 𝑁 𝑇 𝑎(𝑥) 𝑑𝑥
𝑥1 𝑑𝑥
𝑥2
𝐾𝑒3 = ∫ 𝑁 𝑇 𝑏(𝑥) 𝑁𝑑𝑥
𝑥1
𝑥2
𝐹𝑒 = ∫ 𝑁 𝑇 𝑐(𝑥) 𝑑𝑥
𝑥1
L’assemblage pourra se faire suivant l’algorithme décrit au paragraphe 3.8
Discrétisation et application des conditions aux limites :
𝑛+1

𝐾𝑈𝑛 = 𝐹 ⟹ ∑ 𝑘𝑖𝑗 𝑢𝑗 = 𝐹𝑖 , 𝑖 = 1, … , 𝑛 + 1 (𝑖)

𝑗=1
𝑑𝑢 ∆𝑢
En faisant l’approximation de dérivée 𝑑𝑥 ≈ ℎ
les conditions aux limites peuvent s’écrire :
,
𝑢2 − 𝑢1
𝛼1 𝑢1 + 𝛽1 = 𝛾1
{ ℎ
𝑢𝑛+1 − 𝑢𝑛
𝛼2 𝑢𝑛+1 + 𝛽2 = 𝛾2
ℎ
Ou encore
(𝛼 ℎ − 𝛽1 )𝑢1 + 𝛽1 𝑢2 = 𝛾1 ℎ
{ 1
(𝛼2 ℎ + 𝛽2 )𝑢𝑛+1 − 𝛽2 𝑢𝑛 = 𝛾2 ℎ
On peut alors examiner les quatre différents cas suivants :
• 1er cas 𝛼1 ℎ − 𝛽1 = 0 et 𝛼2 ℎ + 𝛽2 = 0
On a nécessairement 𝛽1 ≠ 0 et 𝛽2 ≠ 0 car (𝛼1 , 𝛽1 ) ≠ (0,0) et (𝛼2 , 𝛽2 ) ≠ (0,0). Ainsi :
𝛾 ℎ 𝛾 ℎ
𝑢2 = 𝛽1 et 𝑢𝑛 = − 𝛽2
1 2
Par substitution dans (i), on a :
𝑛+1
𝛾1 ℎ 𝛾2 ℎ
∑ 𝑘𝑖𝑗 𝑢𝑗 = 𝐹𝑖 − 𝑢2 𝑘𝑖2 − 𝑢𝑛 𝑘𝑖𝑛 = 𝐹𝑖 − 𝑘𝑖2 + 𝑘
𝛽1 𝛽2 𝑖𝑛
𝑗=1
𝑗≠2
𝑗≠𝑛
On pourra donc supprimer les colonnes 2 et n, ainsi que les lignes 2 et n de la
matrice de rigidité globale K. Le système résultant se présente comme suit :
𝑘11 𝑘13 ⋯ 𝑘1 𝑛−1 𝑘1 𝑛+1 𝑢1 𝐹1 𝐾12 𝐾1𝑛
𝑘31 𝑘33 ⋯ 𝑘3 𝑛−1 𝑘3 𝑛+1 𝑢3 𝐹3 𝐾32 𝐾3 𝑛
⋮ ⋮ ⋮ ⋮ = − 𝑢2 − 𝑢𝑛
𝑘𝑛−1 1 𝑘𝑛−1 3 ⋯ 𝑘𝑛−1𝑛−1 𝑘𝑛−1 𝑛+1 𝑢𝑛−1 𝐹𝑛−1 𝐾𝑛−1 2 𝐾𝑛−1 𝑛
[𝑘𝑛+1 1 𝑘𝑛+1 3 ⋯ 𝑘𝑛+1𝑛−1 𝑘𝑛+1 𝑛+1 ] [𝑢𝑛+1 ] [𝐹𝑛+1 ] [𝐾𝑛+1 2 ] [𝐾𝑛+1 𝑛 ]
Après la résolution de ce système d’équation linéaire, les composantes 𝑢2 et 𝑢𝑛
du vecteur u, dont les expressions sont connues, pourront être calculées et le
vecteur u pourra être reconstitué.
• 2ème cas : 𝛼1 ℎ − 𝛽1 = 0 et 𝛼2 ℎ + 𝛽2 ≠ 0
On peut écrire alors
𝛾 ℎ 𝛽2 𝛾2 ℎ
𝑢2 = 𝛽1 et 𝑢𝑛+1 = 𝛼 ℎ+𝛽 𝑢𝑛 + 𝛼 ℎ+𝛽
1 2 2 2 2
Et par substitution dans (i), on obtient
𝑛
𝛽2 𝛾2 ℎ
∑ 𝑘𝑖𝑗 𝑢𝑗 + (𝑘𝑖 𝑛 + 𝑘𝑖 𝑛+1 ) 𝑢𝑛 = 𝐹𝑖 − 𝑢2 𝑘𝑖2 − 𝑘
𝛼2 ℎ + 𝛽2 𝛼2 ℎ + 𝛽2 𝑖 𝑛+1
𝑗=1
𝑗≠2
𝑗≠𝑛
On pourra donc supprimer les colonnes 2 et n+1, ainsi que les lignes 2 et n+1 de
la matrice de rigidité globale K. Le système résultant se présente comme suit :

44
 ⋯ ′
𝑘11 𝑘13 𝑘1 𝑛−1 𝑘1𝑛 𝑢1 𝐹1 𝐾12 𝐾1𝑛+1
𝑘31 𝑘33 ⋯ ′
𝑘3 𝑛−1 𝑘3𝑛 𝑢3 𝐹3 𝐾32 𝛾2 ℎ 𝐾3 𝑛+1
⋮ ⋮ ⋮ ⋮ = − 𝑢2 −
⋯ 𝑘𝑛−1𝑛−1 ′ 𝑢𝑛−1 𝐹𝑛−1 𝐾𝑛−1 2 𝛼2 ℎ + 𝛽2 𝐾
𝑘𝑛−11 𝑘𝑛−1 3 𝑘𝑛−1𝑛 𝑛−1 𝑛
[ 𝑘𝑛1 𝑘𝑛 3 ⋯ 𝑘𝑛 𝑛−1 ′
𝑘𝑛 𝑛 ] [ 𝑢 𝑛 ] [ 𝐹𝑛 ] [ 𝐾𝑛 2 ] [𝐾 𝑛 𝑛+1 ]
𝛽2
Où 𝑘𝑖′ 𝑛 = 𝑘𝑖 𝑛 + 𝛼 ℎ+𝛽 𝑘𝑖 𝑛+1 avec 𝑖 = 1,3,4, … , 𝑛
2 2
Après la résolution de ce système d’équation linéaire, les composantes 𝑢2 et 𝑢𝑛+1 du vecteur u,
dont les expressions sont connues (𝑢2 est connue directement, mais 𝑢𝑛+1 est fonction de 𝑢𝑛 ),
pourront être calculées et le vecteur u pourra être reconstitué.
• 3ème cas : 𝛼1 ℎ − 𝛽1 ≠ 0 et 𝛼2 ℎ + 𝛽2 = 0
On a
(𝛼 ℎ − 𝛽1 )𝑢1 + 𝛽1 𝑢2 = 𝛾1 ℎ
{ 1
(𝛼2 ℎ + 𝛽2 )𝑢𝑛+1 − 𝛽2 𝑢𝑛 = 𝛾2 ℎ
𝛽1 𝛾1 ℎ 𝛾 ℎ
𝑢1 = − 𝛼 ℎ−𝛽 𝑢2 + 𝛼 ℎ−𝛽 et 𝑢𝑛 = − 𝛽2
1 1 1 1 2

𝑛+1
𝛽1 𝛾1 ℎ
∑ 𝑘𝑖𝑗 𝑢𝑗 + (𝑘𝑖 2 − 𝑘𝑖 1 ) 𝑢2 = 𝐹𝑖 − 𝑘 − 𝑢𝑛 𝑘𝑖 𝑛
𝛼1 ℎ − 𝛽1 𝛼1 ℎ − 𝛽1 𝑖 1
𝑗=3
𝑗≠𝑛

Et le système à résoudre est le suivant :

′ ⋯ 𝑘2 𝑛−1
𝑘22 𝑘23 𝑘2 𝑛+1 𝑢2 𝐹2 𝐾21 𝐾2 𝑛
′ ⋯ 𝑘3 𝑛−1 𝑘3 𝑛+1
𝑘32 𝑘33 𝑢3 𝐹3 𝛾1 ℎ 𝐾31 𝐾3 𝑛
⋮ ⋮ ⋮ ⋮ = − − 𝑢 𝑛
𝛼1 ℎ − 𝛽1 𝐾
′
𝑘𝑛−12 𝑘𝑛−1 3 ⋯ 𝑘𝑛−1𝑛−1 𝑘𝑛−1 𝑛+1 𝑢𝑛−1 𝐹𝑛−1 𝑛−1 1 𝐾𝑛−1 𝑛
′
[𝑘𝑛+12 𝑘𝑛+1 3 ⋯ 𝑘𝑛+1 𝑛−1 𝑘𝑛+1 𝑛+1 ] 𝑛+1 [𝑢 ] [𝐹𝑛+1 ] [𝐾𝑛+1 1 ] [𝐾𝑛+1 𝑛 ]
′ 𝛽1
Où 𝑘𝑖 2 = 𝑘𝑖 2 − 𝑘𝑖 1 avec 𝑖 = 2,3, … , 𝑛 − 1, 𝑛 + 1
𝛼1 ℎ−𝛽1
Les composantes 𝑢1 et 𝑢𝑛 du vecteur u, dont les expressions sont connues, pourront être calculées
et le vecteur u pourra être reconstitué.
• 4ème cas : 𝛼1 ℎ − 𝛽1 ≠ 0 et 𝛼2 ℎ + 𝛽2 ≠ 0

𝛽1 𝛾1 ℎ 𝛽2 𝛾2 ℎ
𝑢1 = − 𝛼 𝑢2 + 𝛼 et 𝑢𝑛+1 = 𝛼 𝑢𝑛 + 𝛼
1 ℎ−𝛽1 1 ℎ−𝛽1 2 ℎ+𝛽2 2 ℎ+𝛽2

𝑛−1
𝛽1 𝛽2
∑ 𝑘𝑖𝑗 𝑢𝑗 + (𝑘𝑖 2 − 𝑘𝑖 1 ) 𝑢2 + (𝑘𝑖 𝑛 + 𝑘 )𝑢
𝛼1 ℎ − 𝛽1 𝛼2 ℎ + 𝛽2 𝑖 𝑛+1 𝑛
𝑗=3
𝛾1 ℎ 𝛾2 ℎ
= 𝐹𝑖 − 𝑘𝑖 1 − 𝑘
𝛼1 ℎ − 𝛽1 𝛼2 ℎ + 𝛽2 𝑖 𝑛+1
Ce qui donne le système suivant
𝑘22′
𝑘23 ⋯ 𝑘2 𝑛−1 𝑘2′ 𝑛 𝑢2
′
𝑘32 𝑘33 ⋯ 𝑘3 𝑛−1 𝑘3′ 𝑛 𝑢3
⋮ ⋮ ⋮ ⋮
′ ′ 𝑢𝑛−1
𝑘𝑛−1 2 𝑘𝑛−1 3 ⋯ 𝑘𝑛−1𝑛−1 𝑘𝑛−1 𝑛
′
[ 𝑘𝑛 2 𝑘𝑛 3 ⋯ 𝑘𝑛 𝑛−1 𝑘𝑛 𝑛 ] 𝑢𝑛 ]′ [
𝐹2 𝐾21 𝐾2 𝑛+1
𝐹3 𝛾1 ℎ 𝐾31 𝛾2 ℎ 𝐾3 𝑛+1
= − −
𝐹𝑛−1 𝛼1 ℎ − 𝛽1 𝐾 𝛼2 ℎ + 𝛽2 𝐾
𝑛−1 1 𝑛−1 𝑛+1
[ 𝐹𝑛 ] [ 𝐾𝑛 1 ] [ 𝐾𝑛 𝑛+1 ]
𝛽1 𝛽2
Où 𝑘𝑖′ 2 = 𝑘𝑖 2 − 𝛼 ℎ−𝛽 𝑘𝑖 1 et 𝑘𝑖′ 𝑛 = 𝑘𝑖 𝑛 + 𝛼 ℎ+𝛽 𝑘𝑖 𝑛+1 avec 𝑖 = 2,3, … , 𝑛
1 1 2 2

45
 L’implémentation de cet algorithme donne le script suivant :

%==========================================================================
% Programme pour la résolution des problèmes aux limites
% la forme d^2u/dx^2+a(x)du/dx+b(x)u=c(x), x appartenant à ]x0,xn[
% avec des conditions de Robin:
% alpha1*u(x0)+beta1*u'(x0)=gama1
% alpha2*u(xn)+beta2*u'(xn)=gama2
% par la méthode des éléments finis sur
% Fonctions de forme: polynôme d'interpolation
% de Lagrange
% Auteur: GBAGUIDI Thomas Brice
%==========================================================================
clc;
clear all;
% Intervalle
x0=1;xn=10;
% Nombre de subdivision de l'intervalle [x0,xn] pour la MEF
nn=10; % Nombre de noeuds d'interpolation
ne=30; % Nombre d'élément
n=(nn-1)*ne;
%pas
h=(xn-x0)/n;
% Les conditions aux limites
alpha1=0;
beta1=1;
gama1=4*exp(-1);
alpha2=0;
beta2=1;
gama2=-221*exp(-10);
% Nombre de subdivision de l'intervalle [x0,xn] pour les intégrales
Nint=1000;
% Subdivision
X=linspace(x0,xn,n+1);
% Degré polynôme
% Les fonctions a(t) b(t) et c(x)
a=@(t)(1);
b=@(t)(2./t);
c=@(t)(2*(2-1./t)*exp(-t));
% La matrice Ke
Ke=zeros(n+1,n+1);
Fe=zeros(n+1,1);
XX=zeros(nn,ne);
for i=1:ne
k1=(i-1)*(nn-1)+1;k2=(i-1)*(nn-1)+nn;
for jj=1:nn
XX(jj,i)=X(k1+jj-1);
end
dNTdN=@(t)(dLagMEF(XX(:,i),t)'*dLagMEF(XX(:,i),t));
NTdN=@(t)(LagMEF(XX(:,i),t)'*dLagMEF(XX(:,i),t));
NTN=@(t)(LagMEF(XX(:,i),t)'*LagMEF(XX(:,i),t));
ke=@(t)(-dNTdN(t)+a(t)*NTdN(t)+b(t)*NTN(t));
fe=@(t)(LagMEF( XX(:,i),t )'*c(t));
intke=integrale(ke,XX(1,i),XX(nn,i),Nint);
intfe=integrale(fe,XX(1,i),XX(nn,i),Nint);
Ke(k1:k2,k1:k2)=Ke(k1:k2,k1:k2)+intke;
Fe(k1:k2,1)=Fe(k1:k2,1)+intfe;
end
K=Ke;
F=Fe;
test1=alpha1*h-beta1;test2=alpha2*h+beta2;

46
if(test1==0)
u2=gama1*h/beta1;
if(test2==0)
unm1=-gama2*h/beta2;
F=F-u2*Ke(:,2)-unm1*Ke(:,n);
K(:,2)=[];
K(:,n-1)=[];
K(2,:)=[];
F(2)=[];
K(n-1,:)=[];
F(n-1)=[];
u=K\F;
u=[u(1);u2;u(2:n-2);unm1;u(n-1)];
else
F=F-u2*Ke(:,2)-gama2*h/test2*Ke(:,n+1);
K(:,n)=K(:,n)+beta2/test2*K(:,n+1);
K(:,2)=[];
K(:,n)=[];
K(2,:)=[];
F(2)=[];
K(n,:)=[];
F(n)=[];
u=K\F;
un=(beta2*u(end)+gama2*h)/test2;
u=[u(1);u2;u(2:n-1);un];
end
else
if(test2==0)
unm1=-gama2*h/beta2;
F=F-gama1*h/test1*Ke(:,1)-unm1*Ke(:,n);
K(:,2)=K(:,2)-beta1/test1*K(:,1);
K(:,1)=[];
K(:,n-1)=[];
K(1,:)=[];
F(1)=[];
K(n-1,:)=[];
F(n-1)=[];
u=K\F;
u1=(-beta1*u(1)+gama1*h)/test1;
u=[u1;u(1:n-2);unm1;u(n-1)];
else
F=F-gama1*h/test1*Ke(:,1)-gama2*h/test2*Ke(:,n+1);
K(:,n)=K(:,n)+beta2/test2*K(:,n+1);
K(:,2)=K(:,2)-beta1/test1*K(:,1);
K(:,1)=[];
K(:,n)=[];
K(1,:)=[];
F(1)=[];
K(n,:)=[];
F(n)=[];
u=K\F;
u1=(-beta1*u(1)+gama1*h)/test1;
un=(beta2*u(n-1)+gama2*h)/test2;
u=[u1;u(1:n-1);un];
end
end
uexact=(3*X.^2-2*X-1).*exp(-X);
plot(X,u,'ro',X,uexact); % Représentation graphique de la solution
legend('Solution MEF','Solution exacte');
grid on ;

47
Ce programme a été exécuté :
1.4 1.4
Solution MEF d2y/dt 2+dy/dt+2/t*y=2*(2-1/t)exp(-t) Solution MEF
d2y/dt 2+dy/dt+2/t*y=2*(2-1/t)exp(-t)
Solution exacte alpha1*y(t0)+beta1*y'(t0)=gama1 Solution exacte
alpha1*y(t0)+beta1*y'(t0)=gama1 1.2
1.2 alpha2*y(tn)+beta2*y'(tn)=gama2 alpha2*y(tn)+beta2*y'(tn)=gama2

° Nombre d'élément n=10 ° Nombre d'élément n=30

1 ° Nombre de noeuds 1 ° Nombre de noeuds
d'interpolation par élément: d'interpolation par élément:
n=5 n=10
0.8 ° alpha1=0 0.8 ° alpha1=0
° beta1=1 ° beta1=1
° gama1=4/e ° gama1=4/e
0.6 ° alpha2=0 0.6 ° alpha2=0
° beta2=1 ° beta2=1
10
° gama2=-221/e ° gama2=-221/e10
0.4 0.4

0.2 0.2

0 0
1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10
-3
x 10
0.08 10
d2y/dt 2+dy/dt+2/t*y=2*(2-1/t)exp(-t) erreure absolue
d2y/dt 2+dy/dt+2/t*y=2*(2-1/t)exp(-t) erreure absolue
0.07 alpha1*y(t0)+beta1*y'(t0)=gama1 alpha1*y(t0)+beta1*y'(t0)=gama1
alpha2*y(tn)+beta2*y'(tn)=gama2 8 alpha2*y(tn)+beta2*y'(tn)=gama2
0.06 ° Nombre d'élément n=10 ° Nombre d'élément n=30
° Nombre de noeuds ° Nombre de noeuds
0.05 d'interpolation par élément: 6 d'interpolation par élément:
n=5 n=10
° alpha1=0 ° alpha1=0
0.04 ° beta1=1 ° beta1=1
° gama1=4/e 4 ° gama1=4/e
0.03 ° alpha2=0 ° alpha2=0
° beta2=1 ° beta2=1
0.02 ° gama2=-221/e10 2 ° gama2=-221/e10

0.01
0
0

-0.01 -2
1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10

0.35
0.12
d2y/dt 2+dy/dt+2/t*y=2*(2-1/t)exp(-t) erreure relative
d2y/dt 2+dy/dt+2/t*y=2*(2-1/t)exp(-t) erreure relative
alpha1*y(t0)+beta1*y'(t0)=gama1
0.3 alpha1*y(t0)+beta1*y'(t0)=gama1
alpha2*y(tn)+beta2*y'(tn)=gama2 alpha2*y(tn)+beta2*y'(tn)=gama2
0.1
° Nombre d'élément n=10 ° Nombre d'élément n=30
0.25 ° Nombre de noeuds ° Nombre de noeuds
d'interpolation par élément: d'interpolation par élément:
0.08
n=5 n=10
0.2 ° alpha1=0 ° alpha1=0
° beta1=1 ° beta1=1
° gama1=4/e 0.06 ° gama1=4/e
0.15 ° alpha2=0 ° alpha2=0
° beta2=1 ° beta2=1
° gama2=-221/e10 ° gama2=-221/e10
0.04
0.1

0.05 0.02

0 0
1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10

48