from IPython.display import display, Latex
from IPython.core.display import HTML
%reset -f
%matplotlib inline
%autosave 300
from matplotlib.pylab import *
from time import clock
from scipy.sparse import spdiags
from scipy.sparse.linalg import spsolve
from scipy.sparse import lil_matrix, csr_matrix, spdiags
Méthode élément fini¶
On présente la mise en oeuvre des éléments finis $\mathbb P_1$ 1D pour le problème de Robin suivant: $$ \begin {align} -u''(x)=f(x), \quad x\in ]0,L[,\\ -u'(0)+ a_0 u(0)=b_0,\\ u'(L)+ a_1 u(L)=b_1. \end {align} $$ où $f$ est donné dans $L^2(]0,L[)$ et les constantes $a_0$ et $a_1$ sont bien-sûr positives, dont au moins une stricte afin d'assurer l'ellipticité de l'opérateur sur l'espace $H^1(]0,L[)$. Il convient de vérifier en effet que toutes les hypothèses du théorème de Lax-Milgram sont vérifiées pour la formulation faible: $$ \left (\mathcal V\right ) \begin {cases} \text{Trouver } u\in H^1(]0,L[)\text{ tel que }\forall v\in H^1(]0,L[)\\ \int_0^L u'(x)v'(x)\,dx+a_0u(0)v(0)+a_1u(L)v(L)=\int_0^L f(x)v(x)\,dx+b_0v(0)+b_1v(L). \end {cases} $$ On notera $a(u,v)$ l'expression de la forme bilinéaire intervenant dans le membre de gauche de $\left (\mathcal V\right )$ et on notera $l(v)$ la forme linéaire intervenant dans le membre de droite de $\left (\mathcal V\right )$.
Formulation faible approchée¶
On cherche à résoudre le problème approché en dimension finie, $$ \left (\mathcal V_h\right ) \begin {cases} \text{Trouver } u_h\in V_h\text{ tel que }\forall v\in V_h\\ \int_0^L u_h'(x)v'(x)\,dx+a_0u_h(0)v(0)+a_1u_h(L)v(L)=\int_0^L f(x)v(x)\,dx+b_0v(0)+b_1v(L), \end {cases} $$ où $V_h\subset H^1(]0,L[)$ est de dimension finie et $$ V_h=\{ v\in \mathcal C^0(]0,L[)\text{ tel que }v_{|[x_i,x_{i+1}] } \in \mathbb P_1,~\forall 0\le i< N-1\}, $$ avec $x_i=ih$ où $h$ est le pas de la subdivision uniforme de $[0,L]$, soit $h=\frac L {N-1}$.
Base de l'espace d'approximation¶
On construit une base de $V_h$ à l'aide des fonctions chapeaux suivantes $\Phi_i$ pour $0\le i<N$, $$ \begin {cases} \Phi_i(x)=\frac {x-x_{i-1}} h \text { si }x_{i-1}\le x\le x_{i},\\ \Phi_i(x)=\frac {-x+x_{i+1}} h \text { si }x_{i}\le x\le x_{i+1},\\ \Phi_i(x)=0\text { sinon }. \end {cases} $$ On vérifie aisément que $\{\Phi_i \}_{0\le i< N}$ forme une famille libre et génératrice de $V_h$: $$ \forall v\in V_h,\quad v=\sum_{i=0}^{N-1}v(x_i)\Phi_i. $$
Système linéaire associé à $\left (\mathcal V_h\right )$¶
On cherche $u_h$ qui se décompose sur la base $\{\Phi_i \}_{0\le i< N}$ de dimension $N$, $$ u_h(x)=\sum_{i=0}^{N-1}u_i\Phi_i(x),\quad \forall x\in [0,L]. $$ On notera que les inconnues $u_i$ (coefficients scalaires dans la base) satisfont $u_i=u_h(x_i)$, ce qui sera bien pratique pour visualiser la solution sans avoir à reconstruire $u_h$.
Plutôt que de vérifier $\left (\mathcal V_h\right )$ pour tout $v\in V_h$, il suffit de vérifier la formulation variationelle pour tout $\Phi_i$, $0\le i<N$. On a ainsi $N$ équations scalaires pour $N$ inconnues. Soit le système matricielle, $$ AU=F,\text{ avec }U=(u_0,\cdot,u_{N-1})^t, ~F=(l(\Phi_0),\cdot,l(\Phi_{N-1}))^t, \text{ et }A_{ij}=a(\Phi_j,\Phi_i). $$
Assemblage de la matrice¶
Dans le cas présent, il est facile d'identifier chacun des coefficients du système en remarquant que la matrice est tridiagonale compte tenu du support des fonctions de base qui s'interesctent seulement avec la fonction de base précédente et suivante. Il suffit alors de calculer $a(\Phi_i,\Phi_i)$, $a(\Phi_i,\Phi_{i-1})$ et $a(\Phi_i,\Phi_{i+1})$ en traitant à part $i=0$ et $i=N-1$.
Néanmoins, dans un soucis de généralisation, nous allons présenter une méthode qui s'appliquera tout aussi simplement en dimension supérieure. Cette méthode consiste à construire une matrice dite élémentaire qui évalue le calcul intégral de la formulation faible sur un élément fini et l'assemblage de la matrice globale $A$ se faisant alors par somme des contributions venant de chaque élément.
On construit la matrice élémentaire $$ el=\left (\begin {matrix} \int_0^h \varphi_1'(x)\varphi_1'(x)\, dx &\int_0^h\varphi_1'(x)\varphi_2'(x)\, dx\\ \int_0^h \varphi_2'(x)\varphi_1'(x)\, dx &\int_0^h\varphi_2'(x)\varphi_2'(x)\, dx,\\ \end {matrix}\right ) $$ où, quitte à translater l'intervalle de $x_k$, $\varphi_1$ est la restriction de $\Phi_{k}$ à $[x_k,x_{k+1}]$ et $\varphi_2$ est la restriction de $\Phi_{k+1}$ à $[x_k,x_{k+1}]$. Sur $[0,h]$, on a $$ \varphi_1(x)=\frac {h-x}{h}, \quad \varphi_2(x)=\frac {x}{h}. $$ Ainsi, $$ el=\frac 1 {h}\left (\begin {matrix} 1 &-1\\ -1&1\\ \end {matrix}\right ). $$ La clé du remplissage de la matrice globale $A$ repose sur la formule $$ \int_0^L \Phi_i'(x)\Phi_j'(x)\,dx=\sum_{k=0}^{N-2}\int_{x_k}^{x_{k+1}} \Phi_i'(x)\Phi_j'(x)\,dx, $$ mais aussi sur le fait que sur l'intervalle $[x_k,x_{k+1}]$ la contribution des termes $\Phi_i$ et $\Phi_j$ est non trivial lorsque ceux-ci sont les fonctions de base attachés au sommet $x_k$ et $x_{k+1}$. A une translation près, ce sont les fonctions $\varphi_1$ et $\varphi_2$.
Construction du second membre¶
La construction du second membre nécessite le calcul de $$ \int_0^L f(x)\Phi_k(x)\,dx+b_0\Phi_k(0)+b_1\Phi_k(L). $$ Les termes provenant des conditions limites inhomogènes de Robin sont triviaux à calculer et n'apparaissent que si $k=0$ ou si $k=N-1$. Le terme intégral est plus compliqué. On se contentera d'une formule de quadrature pour l'approcher. Une formule des trapèzes basée sur la subdivision $(x_i)_{i=0}^{n-1}$ donne, $$ \int_0^L f(x)\Phi_k(x)\,dx\sim h\left (\frac {f(x_0)\Phi_k(x_0)} 2+\sum_{l=1}^{N-1}f(x_l)\Phi_k(x_l)+\frac {f(x_{N-1})\Phi_k(x_{N-1})} 2\right ). $$ Mais comme $\Phi_k(x_l)=\delta_{kl}$, $$ \begin {cases} F_0=h\frac {f(x_0)} 2 +b_0,\\ F_k=hf(x_k),~0<k<N-1,\\ F_{N-1}=h\frac {f(x_{N-1})} 2 +b_1. \end {cases} $$
N=71
L=7
h=L/(N-1)
#remplissage de la matrice
#le format LIL (list de list) permet de remplir
#les coefficients d'une matrice creuse en accedant à sa position
#ligne colonne
#A est ici initialisée à la matrice nulle de la taille considérée
A =lil_matrix((N, N))
#construction de la matrice élémentaire du terme integral
el=array([[1,-1],[-1,1]])/h
#Parcours des éléments
for k in range(N-1):
S=array([k,k+1])# ce tableau contient les numéros des 2 sommets
for i in range(2):
for j in range(2):
A[S[i],S[j]]=A[S[i],S[j]]+el[i,j]
# la matrice A ici construite ne contient que la contribution
#intégrale, soit celle qu'on aurait eu avec de CL de Neuman homogène
#
#On va rajouter les termes provenant des CL de Robin
a0=30
a1=20
A[0,0]=A[0,0]+a0
A[N-1,N-1]=A[N-1,N-1]+a1
A.toarray()
# on remarque, miracle, que la matrice obtenue par EF P1 est la même
# que celle obtenue par différence fini à un facteur h près
#construction du second membre
b0=200
b1=400
x=linspace(0,L,N)
fx=2*sin(x)#on a choisi f(x)=2sin(x)
F=h*fx
F[0]=F[0]*0.5+b0
F[N-1]=F[N-1]*0.5+b1
#conversion en csr (compressed storage row) pour resolution system
A=A.tocsr()
#resolution systeme creux
U=spsolve(A,F)
plot(x,U)
#exemple avec des points non équirepartis
N=60
L=1
x=linspace(0,L,N)
x=(linspace(0,sqrt(L),N))**2 # beaucoup de point proche de 0 moins proche de 1
h=x[1:N]-x[0:N-1]
#remplissage de la matrice
#le format LIL (list de list) permet de remplir
#les coefficients d'une matrice creuse en accedant à sa position
#ligne colonne
#A est ici initialisée à la matrice nulle de la taille considérée
A =lil_matrix((N, N))
#construction de la matrice élémentaire du terme integral
el=array([[1,-1],[-1,1]])
#Parcours des éléments
for k in range(N-1):
S=array([k,k+1])# ce tableau contient les numéros des 2 sommets
elh=el/h[k]
for i in range(2):
for j in range(2):
A[S[i],S[j]]=A[S[i],S[j]]+elh[i,j]
# la matrice A ici construite ne contient que la contribution
#intégrale, soit celle qu'on aurait eu avec de CL de Neuman homogène
#
#On va rajouter les termes provenant des CL de Robin
a0=200
a1=200
A[0,0]=A[0,0]+a0
A[N-1,N-1]=A[N-1,N-1]+a1
#construction du second membre
b0=200
b1=200
fx=1+0*x#2*sin(x)#on a choisi f(x)=2sin(x)
#Attention à la formule de quadrature qui ne se factorise plus par h constant
F=zeros(N)
for k in range(N-1):
F[k]+=0.5*fx[k]*h[k]
F[k+1]+=0.5*fx[k+1]*h[k]
# robin
F[0]=F[0]+b0
F[N-1]=F[N-1]+b1
#conversion en csr (compressed storage row) pour resolution system
A=A.tocsr()
#resolution systeme creux
U=spsolve(A,F)
plot(x,U,'+')
#Les paramètres du problème de Robin
L=4.5
a0=1
a1=10
# les fonctions choisies pour la validation
u=lambda x: sin(x**2)*sin(x-1)
uprime=lambda x: 2*x*sin(x - 1)*cos(x**2) + sin(x**2)*cos(x - 1)
f=lambda x: 4*x**2*sin(x**2)*sin(x - 1) - 4*x*cos(x**2)*cos(x - 1) + sin(x**2)*sin(x - 1) - 2*sin(x - 1)*cos(x**2)
# on construit b0 et b1 pour satisfaire les conditions limites
# ainsi la fonction u est solution exacte du pb de Robin proposé
b0=-uprime(0)+a0*u(0)
b1=uprime(L)+a1*u(L)
#
# La fonction pour calculer l'erreur
#
def erreur_Robin_P1(N,u,f,p):
# N est le nombre de points de discrétisation
# u et f les fonctions passées en argument, on pourra en changer
# en definissant d'autres fonctions
# p est l'entier pour la norme l^p
h=L/(N-1)
#remplissage de la matrice
#le format LIL (list de list) permet de remplir
#les coefficients d'une matrice creuse en accedant à sa position
#ligne colonne
#A est ici initialisée à la matrice nulle de la taille considérée
A =lil_matrix((N, N))
#construction de la matrice élémentaire du terme integral
el=array([[1,-1],[-1,1]])/h
#Parcours des éléments
for k in range(N-1):
S=array([k,k+1])# ce tableau contient les numéros des 2 sommets
for i in range(2):
for j in range(2):
A[S[i],S[j]]=A[S[i],S[j]]+el[i,j]
#On va rajouter les termes provenant des CL de Robin
A[0,0]=A[0,0]+a0
A[N-1,N-1]=A[N-1,N-1]+a1
#construction du second membre
x=linspace(0,L,N)
F=h*f(x)
F[0]=F[0]*0.5+b0
F[N-1]=F[N-1]*0.5+b1
#conversion en csr (compressed storage row) pour resolution system
A=A.tocsr()
#resolution systeme creux
U=spsolve(A,F)
#le vecteur erreur err
err=abs(U-u(x))
normerr=sum([ei**p for ei in err])
normerr=normerr**(1/p)
#erreur relative
normeU=sum([abs(ui)**p for ui in U])
normeU=normeU**(1/p)
#return max(err)/max(abs(U))
return normerr/normeU
#on evalue l'erreur
print(erreur_Robin_P1(40,u,f,2))
errtab=[]
Ntab=[]
for N in range(20,400,30):
Ntab.append(N)
errtab.append(erreur_Robin_P1(N,u,f,2))
plot(log(Ntab),log(errtab),log(Ntab),-2*log(Ntab))
# on obtient l'ordre 2, c'est mieux que ce que le lemme de Céa garantit
# on parle de superconvergence apportée par la discrétisation régulière.
Validation Dirichlet¶
Par un choix pertient de $a_0$, $a_1$, $b_0$ et $b_1$, proposer avec le même code une validation pour un problème de Dirichlet inhomogène.
Il suffit de choisir $b_0=\beta_0 a_0$, prendre $a_0>>1$ et on obtient formellement $u(0)\sim \beta_0$. On procède de même en $x=L$.
On pose $a_0=a_1=\frac 1 {\epsilon}$ avec $\epsilon>0$, et on note $u_\epsilon$ la solution de notre problème elliptique paramétré par $\epsilon>0$. La preuve rigoureuse de $$ \lim_{\epsilon\to 0^+} u_\epsilon=u, \text{ dans } H^1(0,L), $$ avec $u$ la solution du problème de Dirichlet inhomogène avec $u(0)=\beta_0$, $u(L)=\beta_1$ pourraît très bien faire l'objet d'une partie l'examen.
beta0=2
beta1=11
a0=1.e8
a1=1.e8
b0=1.e8*beta0
b1=1.e8*beta1
N=300
L=7
h=L/(N-1)
A =lil_matrix((N, N))
#construction de la matrice élémentaire du terme integral
el=array([[1,-1],[-1,1]])/h
#Parcours des éléments
for k in range(N-1):
S=array([k,k+1])# ce tableau contient les numéros des 2 sommets
for i in range(2):
for j in range(2):
A[S[i],S[j]]=A[S[i],S[j]]+el[i,j]
A[0,0]=A[0,0]+a0
A[N-1,N-1]=A[N-1,N-1]+a1
#construction du second membre
x=linspace(0,L,N)
fx=2*sin(x)#on a choisi f(x)=2sin(x)
F=h*fx
F[0]=F[0]*0.5+b0
F[N-1]=F[N-1]*0.5+b1
#conversion en csr (compressed storage row) pour resolution system
A=A.tocsr()
#resolution systeme creux
U=spsolve(A,F)
plot(linspace(0,L,N),U)
#on observe le CL beta0 et beta1
# on valide avec les parametres
L=3
beta0=u(0)
beta1=u(L)
a0=1.e8
a1=1.e8
b0=1.e8*beta0
b1=1.e8*beta1
Ntab=[]
errtab=[]
for N in range(20,1000,70):
Ntab.append(N)
errtab.append(erreur_Robin_P1(N,u,f,2))
plot(log(Ntab),log(errtab),log(Ntab),-2*log(Ntab))
Elément fini $\mathbb P_2$¶
Proposer la construction de l'approximation $\mathbb P_2$ conforme de Lagrange pour ce même problème de Robin inhomogène. Vérifier l'ordre de la méthode obtenue. On choisira alors une approximation par formule de Quadrature de Simpson pour le calcul du second membre.
L'élément fini $\mathbb P_2$ conforme est défini à l'aide de trois degrés de liberté, les évaluations aux extrémités de l'élément ainsi qu'au milieu. Les degrés de liberté évalués aux extrémités garantissent le caractère $\mathcal C^0$ de la solution de la méthode élément fini (cf cours).
Formulation approchée¶
On cherche à résoudre le problème approché en dimension finie, $$ \left (\mathcal V_h\right ) \begin {cases} \text{Trouver } u_h\in V_h\text{ tel que }\forall v\in V_h\\ \int_0^L u_h'(x)v'(x)\,dx+a_0u_h(0)v(0)+a_1u_h(L)v(L)=\int_0^L f(x)v(x)\,dx+b_0v(0)+b_1v(L), \end {cases} $$ où $V_h\subset H^1(]0,L[)$ est de dimension finie et $$ V_h=\{ v\in \mathcal C^0(]0,L[)\text{ tel que }v_{|[x_i,x_{i+1}] } \in \mathbb P_2,~\forall 0\le i< N-1\}, $$ avec $x_i=ih$ où $h$ est le pas de la subdivision uniforme de $[0,L]$, soit $h=\frac L {N-1}$.
Base de l'espace d'approximation¶
On construit une base de $V_h$ à l'aide des fonctions suivantes $\Phi_i$ pour $0\le i\le 2N-2$ zt on introduit les points $z_i=i \frac h 2$ pour $0\le i\le 2N-2$. On a $z_{2k}=x_k$ bien-sûr. $$ \begin {cases} {\Phi_{2i}}_{[x_{i-1},x_{i}]}\in \mathbb P_2 \text { interpole } (x_{i-1},0)~(z_{2i-1},0)~(x_{i},1)\text { si }x_{i-1}\le x\le x_{i},\\ (x_{i-1}\le x\le x_{i},\\ {\Phi_{2i}}_{[x_i,x_{i+1}]}\in \mathbb P_2 \text { interpole } (x_{i},1)~(z_{2i+1},0)~(x_{i+1},0)\text { si }x_{i}\le x\le x_{i+1},\\ \Phi_{2i}(x)=0\text { sinon }; {\Phi_{2i+1}}_{[x_{i},x_{i+1}]}\in \mathbb P_2 \text { interpole } (x_{i},0)~(z_{2i+1},1)~(x_{i+1},0)\text { si }x_{i}\le x\le x_{i+1},\\ \Phi_{2i+1}(x)=0\text { sinon }. \end {cases} $$ On vérifie aisément que $\{\Phi_i \}_{0\le i< 2N-1}$ forme une famille libre et génératrice de $V_h$: $$ \forall v\in V_h,\quad v=\sum_{i=0}^{2N-2}v(z_i)\Phi_i. $$
Système linéaire associé à $\left (\mathcal V_h\right )$¶
On cherche $u_h$ qui se décompose sur la base $\{\Phi_i \}_{0\le i< 2N-1}$ de dimension $2N-1$, $$ u_h(x)=\sum_{i=0}^{2N-1}u_i\Phi_i(x),\quad \forall x\in [0,L]. $$ On notera que les inconnues $u_i$ (coefficients scalaires dans la base) satisfont $u_i=u_h(z_i)$, ce qui sera bien pratique pour visualiser la solution sans avoir à reconstruire $u_h$.
Plutôt que de vérifier $\left (\mathcal V_h\right )$ pour tout $v\in V_h$, il suffit de vérifier la formulation variationelle pour tout $\Phi_i$, $0\le i<2N-1$. On a ainsi $2N-1$ équations scalaires pour $2N-1$ inconnues. Soit le système matricielle, $$ AU=F,\text{ avec }U=(u_0,\cdot,u_{2N-2})^t, ~F=(l(\Phi_0),\cdot,l(\Phi_{2N-2}))^t, \text{ et }A_{ij}=a(\Phi_j,\Phi_i). $$
Assemblage de la matrice¶
On doit construit la matrice élémentaire et la matrice globale sera construite par le procédé habituel de remplissage par parcours des éléments. $$ el=\left (\begin {matrix} \int_0^h \varphi_1'(x)\varphi_1'(x)\, dx &\int_0^h\varphi_1'(x)\varphi_2'(x)\, dx&\int_0^h\varphi_1'(x)\varphi_3'(x)\, dx\\ \int_0^h \varphi_2'(x)\varphi_1'(x)\, dx &\int_0^h\varphi_2'(x)\varphi_2'(x)\, dx,&\int_0^h\varphi_2'(x)\varphi_3'(x)\, dx\\ \int_0^h \varphi_3'(x)\varphi_1'(x)\, dx &\int_0^h\varphi_3'(x)\varphi_2'(x)\, dx,&\int_0^h\varphi_3'(x)\varphi_3'(x)\, dx\\ \end {matrix}\right ) $$ où, quitte à translater l'intervalle de $x_k$, $\varphi_1$ est la restriction de $\Phi_{2k}$ à $[x_k,x_{k+1}]$, $\varphi_2$ est la restriction de $\Phi_{2k+1}$ à $[x_k,x_{k+1}]$ et $\varphi_3$ est la restriction de $\Phi_{2k+2}$ à $[x_k,x_{k+1}]$. Sur $[0,h]$, on va construire $(\varphi_i)_{1\le i \le 3}$ à l'aide de Sympy et calculer les intégrales également avec cette librairie. LA fonction interpolate construit le polynome de Lagrange sur les points mentionnés.
import sympy as sy
xs,hs=sy.symbols('xs hs')
phi1=sy.interpolate([(0,1),(hs/2,0),(hs,0)], xs)
phi2=sy.interpolate([(0,0),(hs/2,1),(hs,0)], xs)
phi3=sy.interpolate([(0,0),(hs/2,0),(hs,1)], xs)
phi1p=sy.diff(phi1,xs)
phi2p=sy.diff(phi2,xs)
phi3p=sy.diff(phi3,xs)
print(sy.integrate(phi1p**2,(xs,0,hs)))
print(sy.integrate(phi1p*phi2p,(xs,0,hs)))
print(sy.integrate(phi1p*phi3p,(xs,0,hs)))
print(sy.integrate(phi2p**2,(xs,0,hs)))
print(sy.integrate(phi2p*phi3p,(xs,0,hs)))
print(sy.integrate(phi3p**2,(xs,0,hs)))
On obtient alors la matrice élémentaire $$ el=\frac 1 h \left (\begin {matrix} \frac 7 3 &-\frac 8 3&\frac 1 3\\ -\frac {8} 3 &\frac {16} 3&-\frac 8 3\\ \frac 1 3 &-\frac 8 3&\frac 7 3\\ \end {matrix}\right ) $$
Construction du second membre¶
La construction du second membre nécessite le calcul de $$ \int_0^L f(x)\Phi_k(x)\,dx+b_0\Phi_k(0)+b_1\Phi_k(L). $$ Les termes provenant des conditions limites inhomogènes de Robin sont triviaux à calculer et n'apparaissent que si $k=0$ ou si $k=2N-1$. Le terme intégral est plus compliqué. On se contentera d'une formule de quadrature pour l'approcher. Une formule de Simpson sur chaque élément sera suffisante. $$ \int_0^L f(x)\Phi_k(x)\,dx\sim h \sum_{l=0}^{N-1}\frac 1 6 f(z_{2l})\Phi_k(z_{2l}) +\frac 2 3 f(z_{2l+1})\Phi_k(z_{2l+1})+\frac 1 6 f(z_{2l+2})\Phi_k(z_{2l+2}). $$ On est alors en mesure d'écrire le code.
N=80
L=7
h=L/(N-1)
A =lil_matrix((2*N-1,2*N-1))
#construction de la matrice élémentaire du terme integral
el=array([[7/3.,-8/3.,1/3.],[-8/3.,16/3.,-8/3.],[1/3.,-8/3.,7/3.]])/h
#Parcours des éléments
for k in range(N-1):
S=array([2*k,2*k+1,2*k+2])# ce tableau contient les numéros des 3 ddl
for i in range(3):
for j in range(3):
A[S[i],S[j]]=A[S[i],S[j]]+el[i,j]
# la matrice A ici construite ne contient que la contribution
#intégrale, soit celle qu'on aurait eu avec de CL de Neuman homogène
#
#On va rajouter les termes provenant des CL de Robin
a0=200
a1=100
A[0,0]+=a0
A[2*N-2,2*N-2]+=a1
#construction du second membre
b0=200
b1=200
z=linspace(0,L,2*N-1)
fx=1+2*sin(z)
#formule de quadrature de Simpson
F=zeros(2*N-1)
for k in range(N-1):
F[2*k]+=fx[2*k]*h/6.
F[2*k+1]+=fx[2*k+1]*h*2/3.
F[2*k+2]+=fx[2*k+2]*h/6.
# robin
F[0]+=b0
F[2*N-2]+=b1
#F=h*ones(2*N-1)
#conversion en csr (compressed storage row) pour resolution system
A=A.tocsr()
#resolution systeme creux
U=spsolve(A,F)
plot(z,U)
On va valider la méthode en reprenant le code en fonction.
L=3
beta0=u(0)
beta1=u(L)
a0=1.e8
a1=1.e8
b0=1.e8*beta0
b1=1.e8*beta1
def erreur_Robin_P2(N,u,f,p):
# N est le nombre de points de discrétisation
# u et f les fonctions passées en argument,
# p est l'entier pour la norme l^p
h=L/(N-1)
A =lil_matrix((2*N-1, 2*N-1))
#construction de la matrice élémentaire du terme integral
el=array([[7/3.,-8/3.,1/3.],[-8/3.,16/3.,-8/3.],[1/3.,-8/3.,7/3.]])/h
#Parcours des éléments
for k in range(N-1):
S=array([2*k,2*k+1,2*k+2])# ce tableau contient les numéros des 3 ddl
for i in range(3):
for j in range(3):
A[S[i],S[j]]+=el[i,j]
A[0,0]+=a0
A[2*N-2,2*N-2]+=a1
#construction du second membre
z=linspace(0,L,2*N-1)
fx=f(z)
F=zeros(2*N-1)
for k in range(N-1):
F[2*k]+=fx[2*k]*h/6.
F[2*k+1]+=fx[2*k+1]*h*2/3.
F[2*k+2]+=fx[2*k+2]*h/6.
# robin
F[0]+=b0
F[2*N-2]+=b1
#conversion en csr (compressed storage row) pour resolution system
A=A.tocsr()
#resolution systeme creux
U=spsolve(A,F)
#le vecteur erreur err
err=abs(U-u(z))
#on présente le graphe des erreurs selon les points de discretisation
# plot(x,err)
# print([ei**p for ei in err],sum([ei**p for ei in err]))
normerr=sum([ei**p for ei in err])
normerr=normerr**(1/p)
#erreur relative
normeU=sum([abs(ui)**p for ui in U])
normeU=normeU**(1/p)
#return max(err)/max(abs(U))
return normerr/normeU
# on valide avec les parametres
Ntab=[]
errtab=[]
for N in range(20,300,30):
Ntab.append(N)
errtab.append(erreur_Robin_P2(N,u,f,8))
plot(log(Ntab),log(errtab),log(Ntab),-4*log(Ntab))
#on observe bien une convergence à l'ordre 4 malgré le décrochement
#pour N grand, dû au fait que l'erreur machine devient prépondérante.
