Différences finies centrées

Dans cette méthode , les dérivées première et seconde sont approximées comme suit :
$$
u''_i = \frac{u'_{i+1}-u'_i}{\Delta t}
\quad; \quad u'_i = \frac{u_{i}-u_{i-1}}{\Delta t}
\quad; \quad u'_{i+1} = \frac{u_{i+1}-u_i}{\Delta t}
$$
d'où
\begin{equation}\label{eq:ddui}
u''_i = \frac{u_{i+1}-2u_i +u_{i-1}}{\Delta t^2}
\end{equation}

Cette équation utilise les valeurs de la fonction $u$ aux pas $i+1$,$i$ et $i-1$ pour évaluer la seconde dérivée
au pas $i$. De même on utilisera les pas $i+1$ et $i-1$ pour évaluer la première dérivée au pas $i$.
C'est pour ça qu'on parle de différences finies centrées.
\begin{equation}\label{eq:dui}
u'_i = \frac{u_{i+1}-u_{i-1}}{2\Delta t}
\end{equation}
On remplace maintenant dans l'équation différentielle du second ordre : $M u'' +C u' + K u = P(t)$ et on obtient
$$
M \left( \frac{u_{i+1}-2u_i +u_{i-1}}{\Delta t^2}\right)
+ C \left(\frac{u_{i+1}-u_{i-1}}{2\Delta t}\right)
+ K u_i = P(t_i)
$$
C'ette équation s'écrit pour $u_{i=1}$ comme suit :
\begin{equation}
\left(\frac{1}{\Delta t^2} M + \frac{1}{2\Delta t} C \right) u_{i+1}
= P_i
- \left( K - \frac{2}{\Delta t^2} M \right) u_{i}
- \left(\frac{1}{\Delta t^2} M - \frac{1}{2\Delta t} C \right) u_{i-1}
\end{equation}

sous forme condensée:
\begin{equation}
K^* u_{i+1} = P_i^*
\end{equation}

La solution donne la solution au temps $t_{i+1}$ en fonction des valeurs des deux pas précédents (au temps $t_i$ et $t_{i-1}$).
De ce fait la méthode des différences finies centrées est une méthode explicite.
Au démarrage des calculs, le pas $i=1$ pose problème car la détermination de $u_1$ nécessite non seulement $u_0$ qui est une condition
initiale (une donnée) mais aussi $u_{-1}$ qu'il fait évaluer au préalable.

L'approximation des dérivées (équations \ref{eq:ddui} et \ref{eq:dui}) donne pour $i=0$ :
\begin{align}
u'_0 &= \frac{u_{1}-u_{-1}}{2\Delta t} \\
u''_0 &= \frac{u_{1}-2u_0 + u_{-1}}{(\Delta t)^2}
\end{align}

A partir de ces deux équations, on obtient :
\begin{equation}\label{eq:u_moins1}
u_{-1} = u_0 +\Delta t u'_0 + \frac{\Delta t^2}{2} u''_0
\end{equation}

$u_0$ et $u'_0$ sont les conditions initiales (des donées), reste à déterminer $u''_0$, pour cela
on utilise l'équation de mouvement au tempts $t=t_0$
$$
M u''_0 + C u'_0 +K u_0 = P_0
$$
de laquelle on tire :
\begin{equation}
u''_0 = (P_0 - C u'_0 - K u_0 )/M
\end{equation}

Script MATLAB

function U = dfc(M,C,K,P,Y0,T)

%----------------------------------------------------------------
% Méthode des différences finis centrées
%----------------------------------------------------------------

% initialisation aux valeurs initiales
U0 = Y0(1);
V0 = Y0(2);

% calcul de A0 et U(-1)
dt = T(2)-T(1);
t = T(1);
m = M(t);
c = C(t);
k = K(t);
p = P(t);

A0 = (p - k*U0 - c*V0)/m;
Um1 = U0 - dt*V0 + 0.5*(dt^2)*A0;


% Calcul de U(2) pour amorcer les itérations
Kstar = m/dt^2 + 0.5*c/dt;
a = m/dt^2 - 0.5*c/dt;
b = k - 2*m/dt^2;
Pstar = p - a * Um1 - b*U0;

U(1) = U0;
U(2) = Pstar/Kstar

n = length(T); % nombre de pas de calcul
for i=2:n-1
t=T(i);
dt = T(i+1)-T(i);
m = M(t);
c = C(t);
k = K(t);
p = P(t);
Kstar = m/dt^2 + 0.5*c/dt;
a = m/dt^2 - 0.5*c/dt;
b = k - 2*m/dt^2;
Pstar = p - a*U(i-1) - b*U(i);
U(i+1) = Pstar/Kstar;
end
return

Test de la fonction

%----------------------------------------------------------------
% test de la fonction dfc.m : comparaison avec ode45 de MATLAB
%----------------------------------------------------------------
% Paramètres de l'équation : m u"+ c u' + k u = p
M = @(t) 1+t*exp(-.2*t);
C = @(t) (.1)*sqrt(t);
K = @(t) 1-((1-t))*exp(-.5*t);
P = @(t) sin(t);
% Interval de temps de calcul et conditions initiales
Tb=[0 50];

% Mise sous forme de deux equations : {Y}' = {F(t,Y)} pour ode45
% y'' = (p - k y - c y')/m
% V1 = y ; V2 = y'
% V1' = V2 ; V2'= (p - k y - c y')/m
F = @(t,V) [V(2); (P(t) - K(t)*V(1) - C(t)*V(2))/M(t)]
[T,U] = ode45(F,Tb,Y0,1e-10);
plot(T,U(:,1))
% Calcul avec la fonction dfc
T = [Tb(1):0.5:Tb(2)];
V = dfc(M,C,K,P,Y0,T);
hold on, plot(T,V,'r o')