Niveau Master

Resolution d'edp avec Scilab

Posté par
hesperados 21-08-15 à 21:49

Bonjour,

Je cherche à écrire un programme Scilab pour résoudre en dimension 1 le problème
-u"(x)=f(x) sur [0,1]
u'(0)=0
u'(1)=0
avec des éléments finis.

Comme je n'ai jamais fait d'EDP (je ne suis qu'une informaticienne :p). Je viens vers vous si vous pouviez m'aider à résoudre ce problème.

Posté par re : Resolution d'edp avec Scilab 23-08-15 à 12:07

bonjour : )

pour ton problème je pense que c'est plutôt :
$\large \boxed{\left\{\begin{matrix}-u''(x) = f(x),  x \in [0,1] \\ u(0) = u(1) = 0\end{matrix}\right.}$
du moment que f est continue (ce qui est le cas des fonctions de la vie réelle) la solution au problème est unique et est de classe C² sur [0, 1]

*** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** ***

je te propose ci-dessous une méthode de résolution,
je vais écrire avec une notation Matlab (avec des vecteurs qui ont des indices qui débutent à 1...) _{Scilab ne doit pas être très différent ?}

pour faire général je vais également transformer légèrement le problème comme suit :
on cherche une solution approchée au problème $\large \boxed{\left\{\begin{matrix}-u''(x) + c(x)u(x) = f(x),  x \in [0,1] \\ u(0) = \alpha, u(1) = \beta\end{matrix}\right.}$
(dans ton cas $\boxed{c(x) = 0, \alpha = \beta = 0}$ )

*** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** ***

donc, pour résoudre le problème, on commence par discrétiser notre intervalle [0,1] avec un pas de $\boxed{h = 1/(N + 1)}$ ,
$N$ étant un entier, on notera le vecteur $\boxed{x = (x_1, ..., x_{N+2})^T}$ où $\boxed{x_i = (i - 1)h, i = 1,...,N+2}$

pour plus de commodités on utilisera également les notations suivants : $\boxed{F = (f_1,...,f_{N+2})^T,  f_i = f(x_i), i = 1,...,N+2}$ ,   $\boxed{C = (c_1,...,c_{N+2})^T,  c_i = c(x_i), i = 1,...,N+2}$

ce qu'on cherche c'est une approximation de $\boxed{U = (u(x_1),...,u(x_{N+2}))^T}$ , on notera $\boxed{V = (v_1,...,v_{N+2})^T}$ la solution approchée de U

ici, avec les conditions de bords, on a bien évidemment $\boxed{v_1 = u(0) = \alpha,  v_{N+2} = u(1) = \beta}$ , reste à déterminer $\boxed{W = (v_2,...,v_{N+1})^T}$

je ne vais pas détailler les calculs (sauf si tu les demandes ? _{pour infos ce ne seront que des développements de Taylor jusqu'à l'ordre 6}), mais le problème approché se résume à la résolution du produit matriciel suivant $\boxed{(P_h) : AW = b}$ avec :

$\large \boxed{A = A_1 + \frac{h^2}{12}A_2 \in \mathbb{R}^{N\times N}}   \boxed{b = \frac{h^2}{12}b_1 + b_2 \in \mathbb{R}^N} \\ \\ \boxed{A_1 = \begin{pmatrix}2 & -1 &  & 0 \\ -1 & \ddots & \ddots & \\ & \ddots & \ddots & -1 \\ 0 &  & -1 & 2 \\ \end{pmatrix} \in \mathbb{R}^{N\times N}}   \boxed{A_2 = \begin{pmatrix}10c_2 & \, \: \! c_3 &  & 0\\ \\ c_2 & 10c_3 & \ddots & \\ \\ & \ddots & \ddots & c_{N+1}\\ \\ 0 &  & c_N & 10c_{N+1} \\ \end{pmatrix} \in \mathbb{R}^{N\times N}} \\ \\ \boxed{b_1 = \frac{h^2}{12}\begin{pmatrix}f_1 + 10f_2 + f_3 \\ \vdots \\ \vdots \\ f_N + 10f_{N+1} + f_{N+2} \\ \end{pmatrix} \in \mathbb{R}^N}   \boxed{b_2 = \begin{pmatrix}(1 - \frac{h^2}{12}c_1)\alpha \\ 0 \\ \vdots \\ 0 \\ (1 - \frac{h^2}{12}c_{N+2})\beta \\ \end{pmatrix} \in \mathbb{R}^N}$

la solution approchée est donnée par : $\large \boxed{V = \begin{pmatrix}\alpha \\ W \\ \beta \\ \end{pmatrix}}$

il s'agit d'une méthode d'ordre 4, le schéma de Numerov

Posté par re : Resolution d'edp avec Scilab 23-08-15 à 12:26

juste une petite correction au niveau de b1 il faut enlever le h^2/12, il s'agit de :

$\large \boxed{b = \frac{h^2}{12}b_1 + b_2 \in \mathbb{R}^N} \\ \\ \boxed{b_1 = \begin{pmatrix}f_1 + 10f_2 + f_3 \\ \vdots \\ \vdots \\ f_N + 10f_{N+1} + f_{N+2}\end{pmatrix} \in \mathbb{R}^N} \boxed{b_2 = \begin{pmatrix}(1 - \frac{h^2}{12}c_1)\alpha \\ 0 \\ \vdots \\ 0 \\ (1 - \frac{h^2}{12}c_{N+2})\beta \\ \end{pmatrix} \in \mathbb{R}^N}$