On considère l'équation différentielle du second ordre suivante :
\begin{equation}
y''(x) = p(x) y'(x) + q(x) y(x) + r(x)
\label{ode_ordre_2}
\end{equation}
avec les conditions aux limites :
$$ y(a) = y_a \quad \text{et} \quad y(b)=y_b $$
La discrétisation des dérivées se fait en utilisant deux pas. Nous avons :
\begin{eqnarray}
y(x_i + h) = y(x_i) + h y'(x_i) + h^2 y''(x_i) \\
y(x_i - h) = y(x_i) - h y'(x_i) + h^2 y''(x_i)
\end{eqnarray}
où : $$ y(x_i) =y_i , \quad y(x_i+h) = y_{i+1} , \quad y(x_i-h) = y_{i-1} $$
La somme des deux expressions donne
\begin{equation}
y''(x_i) = y''_i = \frac{y_{i+1} -2 y_i +y_{i-1}}{h^2}
\end{equation}
La dérivée d'ordre 1 est évaluée entre $x_i+h$ et $x_i-h$
\begin{equation}
y'_i = \frac{y_{i+1}-y_{i-1}}{2h}
\end{equation}
En remplace maintenant dans l'équation (\ref{ode_ordre_2}) :
$$
y_{i+1} - 2 y_i +y_{i-1}
=
\frac{h}{2} p(x_i) \left( y_{i+1} - y_{i-1} \right) + h^2 q(x_i) y_i +h^2 r(x_i)
$$
Soit :
\begin{equation}
- \left( 1 + \tfrac{1}{2} h p_i \right) y_{i-1}
+ \left( 2 + h^2 q_i \right) y_i
+ \left( \tfrac{1}{2} h p_i -1 \right) y_{i+1}
= -h^2 r_i
\end{equation}
avec
$$
p_i = p(x_i) , \quad q_i=q(x_i) \quad \text{et}\quad r_i=r(x_i)
$$
Pour $i=1$ :
$$
- \left( 1 + \tfrac{1}{2} h p_1 \right) y_{0}
+ \left( 2 + h^2 q_1 \right) y_1
+ \left( \tfrac{1}{2} h p_1 - 1 \right) y_{2}
= -h^2 r_1
$$
où $y_0 = y(x_0) = y(a) = y_a$
est une donnée. L'équation devient :
\begin{equation}
\left( 2 + h^2 q_1 \right) y_1
+ \left( \tfrac{1}{2} h p_1 - 1 \right) y_{2}
= -h^2 r_1+ \left( 1 + \tfrac{1}{2} h p_1 \right) y_{0}
\end{equation}
De même pour $i=n$, tenant compte de $x_{n+1} = b$ ; $y_{n+1} = y_b$ comme valeur aux limites :
\begin{equation}
- \left( 1 + \tfrac{1}{2} h p_n \right) y_{n-1}
+ \left( 2 + h^2 q_1 \right) y_n
= -h^2 r_n - \left( \tfrac{1}{2} h p_n - 1 \right) y_b
\end{equation}
Sous forme matricielle :
\begin{equation}
\left[\mathbf{A}\right] \left\{\mathbf{y}\right\} = \left\{\mathbf{B}\right\}
\end{equation}
avec
$$
\left[\mathbf{A}\right] =
\begin{bmatrix}
2 + h^2 q_1 & {-1 +\tfrac{1}{2} h p_1} & \qquad\ 0 \qquad \cdots &\cdots\ 0 \qquad \\
{-1-\tfrac{1}{2}h p_2} & {2+ h^2 q_2 } & -1+\tfrac{1}{2} h p_2 \ \ddots & \vdots \\
\qquad\quad 0 \qquad \ddots & \qquad\ddots & \qquad\ddots & \vdots \\
\qquad \vdots \qquad\ & \ddots \qquad\qquad & \ddots \qquad\qquad\qquad \ddots & \qquad 0 \qquad\\
\qquad \vdots \qquad\ &{-1-\tfrac{1}{2}h p_{n-1}} & 2+ h^2 q_{n-1} &{-1+\tfrac{1}{2} h p_{n-1}} \\
\quad 0 \ \cdots & \cdots 0 &-1-\tfrac{1}{2} h p_n & 2+ h^2 q_n
\end{bmatrix}
$$
$$
\left\{\mathbf{y}\right\} =
\begin{Bmatrix}
y_1 \\ y_2 \\ \vdots \\ y_{n-1} \\ y_n
\end{Bmatrix}
\qquad; \qquad
\left\{\mathbf{B}\right\}
\begin{Bmatrix}
- h^2 r_1 +\left(1+\tfrac{1}{2} h p_1\right)y_a \\
- h^2 r_2 \\
\vdots \\
- h^2 r_{n-1} \\
- h^2 r_n + \left(1-\tfrac{1}{2} h p_n\right)y_b
\end{Bmatrix}
$$
La matrice $\mathbf{A}$ du système d'équations algébriques précédent est une matrice tridiagonale. La solution peut être obtenue en réduisant les équations l'une après l'autre.
On note par $L$, $D$ et $S$ les diagonales inférieure, centrale et supérieure de la matrice $\mathbf{A}$, respectivement. On a
$$
\begin{aligned}
L_i &= A_{i+1,i} &;\ &i=1,n-1\\
D_i &= A_{i,i} &;\ &i=1,n \\
S_i &= A_{i-1,i} &;\ &i=2,n
\end{aligned}
$$
Les $n$ equations s'écrivent comme suit :
\begin{eqnarray*}
D_1 y_1 + S_2 y_2 &=& B_1 \\
\cdots \cdots &\cdots& \cdots \\
L_{i-1} y_{i-1} + D_i y_i + S_{i+1} y_{i+1} &=& B_i \\
\cdots \cdots &\cdots& \cdots \\
L_{n-1} y_{n-1} + D_n y_n &=& B_n
\end{eqnarray*}
La première équation permet de déduire $y_1$ du système :
$$ y_1 = B_1/D_1 - (S_2/D_1) y_2$$
En remplaçant dans la deuxième équation :
$$L_1 \left( B_1/D_1 - (S_2/D_1) y_2 \right) + D_2 y_2 + S_3 y_3 = B_2$$
elle se met sous la forme suivante :
$$ D_2^* y_2 + S_3 y_3 = B_2^*$$
avec
$$
D_2^* = D_2 - \frac{L_1}{D_1}S_2 \quad;\quad
B_2^* = B_2 - \frac{L_1 B_1}{D_1}
$$
Cette équation permet à son tour de déduire $y_2$ :
$ y_2 = B_2^*/D_2^* - (S_3/D_2^*) y_3$
et de le remplacer dans la troisième équation qui se met sous la même forme que celle d'avant :
$$
D_3^* y_3 + S_4 y_4 = B_3^* \quad;\quad
D_3^* = D_3 - \frac{L_2}{D_2^*}S_3\ , \
B_3^* = B_3 - L_2 \frac{B_2^*}{D_2^*}
$$
De cette manière on écrit la $i^\text{ème}$ équation :
$$ D_i^* y_i + S_{i+1} y_{i+1} = B_i^*$$
d'où on calcul $y_i$ par :
\begin{equation}
y_i = \frac{B_i^*}{D_i^* } - \frac{S_{i+1}}{D_i^*} y_{i+1}
\label{eq:tridiag:yi}
\end{equation}
avec
\begin{equation}
D_i^* = D_i - \frac{L_{i-1}}{D_{i-1}^*} S_i \quad;\quad
B_i^* = B_i - L_{i-1} \frac{B_{i-1}^*}{D_{i-1}^*}
\label{eq:tridiag:DiBi}
\end{equation}
La dernière équation se réduit à :
$$
D_n^* y_n = B_n^* \quad;\quad
D_n^* = D_n - \frac{L_{n-1}}{D_{n-1}^*}S_n\ , \
B_n^* = B_n - L_{n-1} \frac{B_{n-1}^*}{D_{n-1}^*}
$$
Ce qui permet de calculer $y_n$, puis $y_{n-1}$, $y_{n-2}$, $y_{n-3}$, ..., $y_2$ et $y_1$ en utilisant l'expression de récurrence (\ref{eq:tridiag:yi}).
Comme exemple, on considère l'équation $y'' = x - y$ avec
les conditions $y(0) = 0$ et $y(5)=6 $.
La solution exacte est $y_e = x + \sin(x)/\sin(5) \approx x \sin(x)$.
Les fonctions de $x$ dans l'équation (\ref{ode_ordre_2}) sont dans ce cas :
$$p(x)=0 \ , \quad q(x)=-1 \ , \quad r(x) = x$$
L'intervalle $[0,5]$ est divisé en $10$ segments de taille $h=0.5$ chacun.
Les coordonnées des points sont $x_0=0, x_1=0.5, x_2=1 \cdots , x_9=5.5, x_{10}=5$
Les composantes de la matrice $\mathbf{A}$ prennent les valeurs suivantes :
\begin{eqnarray*}
A_{ii}&=&2+0.5^2 \times (-1) = 1.75 \\
A_{ii+1}&=& -1 + 0.5\times0.5 \times 0 = -1 \\
A_{i+1i}&=& -1 - 0.5\times0.5 \times 0 = -1
\end{eqnarray*}
Les composantes du vecteur $\mathbf{B}$ prennent les valeurs :
\begin{eqnarray*}
B_1 &=& -0.5^2 \times 0.5 +(1+0.5\times 0)\times 0 = -0.125\\
B_2 &=& -0.5^2 \times 1 = -0.25, \quad \\
B_3 &=& -0.5^2 \times 1.5 = -0.375 \\
\cdots&\cdots& \cdots \cdots \\
B_9 &=& -0.5^2 \times 4.5 +(1-0.5\times 0)\times 6 = 4.875
\end{eqnarray*}
Le système algébrique est :
$$
\begin{bmatrix}
1.75 & -1 & 0 & \cdots & \cdots & & & & 0 \\
-1 & 1.75 & -1 &0 &&&&&\vdots\\
0 & -1 & 1.75& -1&0 &&&&\vdots\\
\vdots & \ddots & -1 & 1.75&-1&0 \\
& & & -1 & 1.75&-1&0 \\
&&& & -1 & 1.75&-1 &0 \\
&&&& & -1 & 1.75&-1 &0 \\
\vdots &&&&& \ddots & -1 & 1.75 &-1\\
0 &\cdots&&&&\cdots& 0& -1 & 1.75 \\
\end{bmatrix}
\begin{Bmatrix}
y_1 \\ y2\\ \vdots \\ \\ \\ \\ \vdots\\ y_8\\ y_9
\end{Bmatrix}
=
\begin{Bmatrix}
-0.125\\
-0.25\\
-0.375\\
-0.5\\
-0.625\\
-0.75\\
-0.875\\
-1\\
4.875
\end{Bmatrix}
$$
b=5;
ya=0;
yb=6;
n = 9;
h = (b-a)/(n+1);
p = @(x) 0*x;
q = @(x) -1+0*x;
r = @(x) x;
x = [1:n]*h;
L = -1-h/2*p(x(2:n));
D = 2+(h^2)*q(x);
S = -1+(h/2)*p(x(1:n-1));
B = -h^2*r(x');
B(1) = B(1)+(1+(h/2)*p(x(1)))*ya;
B(n) = B(n)+(1-(h/2)*p(x(n)))*yb;
% Solution du système tridiagonal
u = 0*D;
S=[0 S];
for i=2:n % equation (\ref{eq:tridiag:DiBi})
D(i) = D(i) - L(i-1)*S(i)/D(i-1);
B(i) = B(i) - L(i-1)*B(i-1)/D(i-1);
end
u(n) = B(n)/D(n);
for i=n-1:-1:1
u(i) = B(i)/D(i) - S(i+1)/D(i) * u(i+1);
end
u = [ya u yb]; % ajout des CALs
x = [a x b];
ue = x+sin(x)/sin(5);
plot(x,ue,'r', x,u,'o')
\begin{equation}
y''(x) = p(x) y'(x) + q(x) y(x) + r(x)
\label{ode_ordre_2}
\end{equation}
avec les conditions aux limites :
$$ y(a) = y_a \quad \text{et} \quad y(b)=y_b $$
La discrétisation des dérivées se fait en utilisant deux pas. Nous avons :
\begin{eqnarray}
y(x_i + h) = y(x_i) + h y'(x_i) + h^2 y''(x_i) \\
y(x_i - h) = y(x_i) - h y'(x_i) + h^2 y''(x_i)
\end{eqnarray}
où : $$ y(x_i) =y_i , \quad y(x_i+h) = y_{i+1} , \quad y(x_i-h) = y_{i-1} $$
La somme des deux expressions donne
\begin{equation}
y''(x_i) = y''_i = \frac{y_{i+1} -2 y_i +y_{i-1}}{h^2}
\end{equation}
La dérivée d'ordre 1 est évaluée entre $x_i+h$ et $x_i-h$
\begin{equation}
y'_i = \frac{y_{i+1}-y_{i-1}}{2h}
\end{equation}
En remplace maintenant dans l'équation (\ref{ode_ordre_2}) :
$$
y_{i+1} - 2 y_i +y_{i-1}
=
\frac{h}{2} p(x_i) \left( y_{i+1} - y_{i-1} \right) + h^2 q(x_i) y_i +h^2 r(x_i)
$$
Soit :
\begin{equation}
- \left( 1 + \tfrac{1}{2} h p_i \right) y_{i-1}
+ \left( 2 + h^2 q_i \right) y_i
+ \left( \tfrac{1}{2} h p_i -1 \right) y_{i+1}
= -h^2 r_i
\end{equation}
avec
$$
p_i = p(x_i) , \quad q_i=q(x_i) \quad \text{et}\quad r_i=r(x_i)
$$
Pour $i=1$ :
$$
- \left( 1 + \tfrac{1}{2} h p_1 \right) y_{0}
+ \left( 2 + h^2 q_1 \right) y_1
+ \left( \tfrac{1}{2} h p_1 - 1 \right) y_{2}
= -h^2 r_1
$$
où $y_0 = y(x_0) = y(a) = y_a$
est une donnée. L'équation devient :
\begin{equation}
\left( 2 + h^2 q_1 \right) y_1
+ \left( \tfrac{1}{2} h p_1 - 1 \right) y_{2}
= -h^2 r_1+ \left( 1 + \tfrac{1}{2} h p_1 \right) y_{0}
\end{equation}
De même pour $i=n$, tenant compte de $x_{n+1} = b$ ; $y_{n+1} = y_b$ comme valeur aux limites :
\begin{equation}
- \left( 1 + \tfrac{1}{2} h p_n \right) y_{n-1}
+ \left( 2 + h^2 q_1 \right) y_n
= -h^2 r_n - \left( \tfrac{1}{2} h p_n - 1 \right) y_b
\end{equation}
Sous forme matricielle :
\begin{equation}
\left[\mathbf{A}\right] \left\{\mathbf{y}\right\} = \left\{\mathbf{B}\right\}
\end{equation}
avec
$$
\left[\mathbf{A}\right] =
\begin{bmatrix}
2 + h^2 q_1 & {-1 +\tfrac{1}{2} h p_1} & \qquad\ 0 \qquad \cdots &\cdots\ 0 \qquad \\
{-1-\tfrac{1}{2}h p_2} & {2+ h^2 q_2 } & -1+\tfrac{1}{2} h p_2 \ \ddots & \vdots \\
\qquad\quad 0 \qquad \ddots & \qquad\ddots & \qquad\ddots & \vdots \\
\qquad \vdots \qquad\ & \ddots \qquad\qquad & \ddots \qquad\qquad\qquad \ddots & \qquad 0 \qquad\\
\qquad \vdots \qquad\ &{-1-\tfrac{1}{2}h p_{n-1}} & 2+ h^2 q_{n-1} &{-1+\tfrac{1}{2} h p_{n-1}} \\
\quad 0 \ \cdots & \cdots 0 &-1-\tfrac{1}{2} h p_n & 2+ h^2 q_n
\end{bmatrix}
$$
$$
\left\{\mathbf{y}\right\} =
\begin{Bmatrix}
y_1 \\ y_2 \\ \vdots \\ y_{n-1} \\ y_n
\end{Bmatrix}
\qquad; \qquad
\left\{\mathbf{B}\right\}
\begin{Bmatrix}
- h^2 r_1 +\left(1+\tfrac{1}{2} h p_1\right)y_a \\
- h^2 r_2 \\
\vdots \\
- h^2 r_{n-1} \\
- h^2 r_n + \left(1-\tfrac{1}{2} h p_n\right)y_b
\end{Bmatrix}
$$
Résolution du système d'équations
On note par $L$, $D$ et $S$ les diagonales inférieure, centrale et supérieure de la matrice $\mathbf{A}$, respectivement. On a
$$
\begin{aligned}
L_i &= A_{i+1,i} &;\ &i=1,n-1\\
D_i &= A_{i,i} &;\ &i=1,n \\
S_i &= A_{i-1,i} &;\ &i=2,n
\end{aligned}
$$
Les $n$ equations s'écrivent comme suit :
\begin{eqnarray*}
D_1 y_1 + S_2 y_2 &=& B_1 \\
\cdots \cdots &\cdots& \cdots \\
L_{i-1} y_{i-1} + D_i y_i + S_{i+1} y_{i+1} &=& B_i \\
\cdots \cdots &\cdots& \cdots \\
L_{n-1} y_{n-1} + D_n y_n &=& B_n
\end{eqnarray*}
La première équation permet de déduire $y_1$ du système :
$$ y_1 = B_1/D_1 - (S_2/D_1) y_2$$
En remplaçant dans la deuxième équation :
$$L_1 \left( B_1/D_1 - (S_2/D_1) y_2 \right) + D_2 y_2 + S_3 y_3 = B_2$$
elle se met sous la forme suivante :
$$ D_2^* y_2 + S_3 y_3 = B_2^*$$
avec
$$
D_2^* = D_2 - \frac{L_1}{D_1}S_2 \quad;\quad
B_2^* = B_2 - \frac{L_1 B_1}{D_1}
$$
Cette équation permet à son tour de déduire $y_2$ :
$ y_2 = B_2^*/D_2^* - (S_3/D_2^*) y_3$
et de le remplacer dans la troisième équation qui se met sous la même forme que celle d'avant :
$$
D_3^* y_3 + S_4 y_4 = B_3^* \quad;\quad
D_3^* = D_3 - \frac{L_2}{D_2^*}S_3\ , \
B_3^* = B_3 - L_2 \frac{B_2^*}{D_2^*}
$$
De cette manière on écrit la $i^\text{ème}$ équation :
$$ D_i^* y_i + S_{i+1} y_{i+1} = B_i^*$$
d'où on calcul $y_i$ par :
\begin{equation}
y_i = \frac{B_i^*}{D_i^* } - \frac{S_{i+1}}{D_i^*} y_{i+1}
\label{eq:tridiag:yi}
\end{equation}
avec
\begin{equation}
D_i^* = D_i - \frac{L_{i-1}}{D_{i-1}^*} S_i \quad;\quad
B_i^* = B_i - L_{i-1} \frac{B_{i-1}^*}{D_{i-1}^*}
\label{eq:tridiag:DiBi}
\end{equation}
La dernière équation se réduit à :
$$
D_n^* y_n = B_n^* \quad;\quad
D_n^* = D_n - \frac{L_{n-1}}{D_{n-1}^*}S_n\ , \
B_n^* = B_n - L_{n-1} \frac{B_{n-1}^*}{D_{n-1}^*}
$$
Ce qui permet de calculer $y_n$, puis $y_{n-1}$, $y_{n-2}$, $y_{n-3}$, ..., $y_2$ et $y_1$ en utilisant l'expression de récurrence (\ref{eq:tridiag:yi}).
Exemple
les conditions $y(0) = 0$ et $y(5)=6 $.
La solution exacte est $y_e = x + \sin(x)/\sin(5) \approx x \sin(x)$.
Les fonctions de $x$ dans l'équation (\ref{ode_ordre_2}) sont dans ce cas :
$$p(x)=0 \ , \quad q(x)=-1 \ , \quad r(x) = x$$
L'intervalle $[0,5]$ est divisé en $10$ segments de taille $h=0.5$ chacun.
Les coordonnées des points sont $x_0=0, x_1=0.5, x_2=1 \cdots , x_9=5.5, x_{10}=5$
Les composantes de la matrice $\mathbf{A}$ prennent les valeurs suivantes :
\begin{eqnarray*}
A_{ii}&=&2+0.5^2 \times (-1) = 1.75 \\
A_{ii+1}&=& -1 + 0.5\times0.5 \times 0 = -1 \\
A_{i+1i}&=& -1 - 0.5\times0.5 \times 0 = -1
\end{eqnarray*}
Les composantes du vecteur $\mathbf{B}$ prennent les valeurs :
\begin{eqnarray*}
B_1 &=& -0.5^2 \times 0.5 +(1+0.5\times 0)\times 0 = -0.125\\
B_2 &=& -0.5^2 \times 1 = -0.25, \quad \\
B_3 &=& -0.5^2 \times 1.5 = -0.375 \\
\cdots&\cdots& \cdots \cdots \\
B_9 &=& -0.5^2 \times 4.5 +(1-0.5\times 0)\times 6 = 4.875
\end{eqnarray*}
Le système algébrique est :
$$
\begin{bmatrix}
1.75 & -1 & 0 & \cdots & \cdots & & & & 0 \\
-1 & 1.75 & -1 &0 &&&&&\vdots\\
0 & -1 & 1.75& -1&0 &&&&\vdots\\
\vdots & \ddots & -1 & 1.75&-1&0 \\
& & & -1 & 1.75&-1&0 \\
&&& & -1 & 1.75&-1 &0 \\
&&&& & -1 & 1.75&-1 &0 \\
\vdots &&&&& \ddots & -1 & 1.75 &-1\\
0 &\cdots&&&&\cdots& 0& -1 & 1.75 \\
\end{bmatrix}
\begin{Bmatrix}
y_1 \\ y2\\ \vdots \\ \\ \\ \\ \vdots\\ y_8\\ y_9
\end{Bmatrix}
=
\begin{Bmatrix}
-0.125\\
-0.25\\
-0.375\\
-0.5\\
-0.625\\
-0.75\\
-0.875\\
-1\\
4.875
\end{Bmatrix}
$$
Script Matlab
b=5;
ya=0;
yb=6;
n = 9;
h = (b-a)/(n+1);
p = @(x) 0*x;
q = @(x) -1+0*x;
r = @(x) x;
x = [1:n]*h;
L = -1-h/2*p(x(2:n));
D = 2+(h^2)*q(x);
S = -1+(h/2)*p(x(1:n-1));
B = -h^2*r(x');
B(1) = B(1)+(1+(h/2)*p(x(1)))*ya;
B(n) = B(n)+(1-(h/2)*p(x(n)))*yb;
% Solution du système tridiagonal
u = 0*D;
S=[0 S];
for i=2:n % equation (\ref{eq:tridiag:DiBi})
D(i) = D(i) - L(i-1)*S(i)/D(i-1);
B(i) = B(i) - L(i-1)*B(i-1)/D(i-1);
end
u(n) = B(n)/D(n);
for i=n-1:-1:1
u(i) = B(i)/D(i) - S(i+1)/D(i) * u(i+1);
end
u = [ya u yb]; % ajout des CALs
x = [a x b];
ue = x+sin(x)/sin(5);
plot(x,ue,'r', x,u,'o')
| $x$ | $y_\text{numérique}$ | $y_\text{exacte}$ |
| 0.5 | -0.0137 | 0.0000 |
| 1.0 | 0.1010 | 0.1225 |
| 1.5 | 0.4404 | 0.4598 |
| 2.0 | 1.0448 | 1.0518 |
| 2.5 | 1.8879 | 1.8759 |
| 3.0 | 2.8841 | 2.8528 |
| 3.5 | 3.9093 | 3.8658 |
| 4.0 | 4.8321 | 4.7892 |
| 4.5 | 5.5469 | 5.5194 |
% Exemple : y"(x) = x-y(x)
% CAL : y(0)=0 et y(5)=6