FEM: singularité de la matrice de rigidité

11

Je résous l'équation différentielle avec les conditions initiales u (0) = u (1) = 0 , u '' (0) = u '' (1) = 0 . Ici \ sigma (x) \ geqslant \ sigma_ {0}> 0 est le paramètre. Sous forme d'opérateur, nous pouvons réécrire l'équation différentielle comme Au = f , où l'opérateur A est défini positif.

(σ2(x)u(x))=f(x),0x1
u(0)=u(1)=0u(0)=u(1)=0σ(x)σ0>0Au=fA

En suivant le schéma FEM, je réduit mon problème à un problème d'optimisation J (u) = (Au, u) - 2 (f, u) \ to \ min_ {u} J'introduis les

J(u)=(Au,u)2(f,u)minu
éléments finis hk(x) comme
vk(x)={1(xxkh)2,x[xk1,xk+1]0,otherwise
pour tout k=1,,n1 , où xk=hk , h=1n . Les éléments finis v0(x) et vn(x) sont introduits de manière similaire.

J'essaie de trouver numériquement le vecteur α tel que u(x)=k=0nαkvk(x) résout le problème d'optimisation. Nous avons

J(u)=i=0nj=0nαiαj(Avi,vj)i=0n2αi(vi,f)=αTVα2αTbminα,
bi=(f,vi) et Vi,j=(Avi,vj) . Après différenciation par rapport à α je reçois
Vα=b,
mais ici la matrice de rigidité V est singulière. Alors que dois-je faire? Peut-être que je dois choisir d'autres éléments finis?
Appliqué
la source
Salut Nimza, avez-vous un problème de test dont vous connaissez la solution exacte? Si oui, essayez d'abord de résoudre VTVα=VTb pour tester si votre base est correcte à l'intérieur du domaine, si tout semble correct, alors c'est peut-être le BC incorrectement posé qui rend la matrice singulière. Mais le BC me semble bien.
Shuhao Cao

Réponses:

13

Par ordre décroissant de probabilité

  1. Base incorrecte. D'après votre description, il semble que vous ayez exactement deux fonctions quadratiques avec support sur chaque élément. Cet espace n'est pas une partition d'unité et n'est pas (dérivées premières continues). Pour discrétiser directement votre problème du quatrième ordre (au lieu de le réduire à un système d'équations du deuxième ordre, par exemple), vous aurez besoin d'une base . Notez que la base devrait être capable de reproduire exactement toutes les fonctions linéaires.C1C1C1

  2. Conditions aux limites insuffisantes. Cela sera manifestement évident si vous calculez et tracez l'espace nul.

  3. Assemblage incorrect. Vérifiez la carte des éléments à l'ordre assemblé pour confirmer qu'elle correspond à ce que vous attendiez, par exemple qu'elle n'inverse pas l'orientation des éléments.

  4. Assemblage local incorrect. Dans 1D, vous pouvez calculer analytiquement à quoi ressemble la matrice de rigidité des éléments (peut-être pour un cas simplifié) et vérifier que le code la reproduit.

Jed Brown
la source
Je vous remercie. 1. Je pense que j'aurai besoin d'une base parce que . Ensuite, si je considère uniquement les fonctions qui satisfont aux conditions aux limites, alors . C2(Au,v)=01σ2(x)u(x)v(x)dxkerA={0}
Appliqué
1
Une base suffit, l'intégrande n'a pas besoin d'être continue. Notez que les conditions aux limites sur les dérivées secondes deviendront une intégrale aux limites. Vous pouvez utiliser une base pour discrétiser directement un problème de quatrième ordre, mais vous devrez intégrer des termes de saut comme avec les méthodes Galerkin discontinues pour les systèmes de premier et de second ordre. Ce n'est pas une mauvaise méthode, mais elle est inutilement compliquée dans 1D car il est si facile de construire des bases avec n'importe quel ordre de continuité (par exemple des splines). Ce document est un exemple de " DG". C1C0C0
Jed Brown
D'accord. J'ai corrigé ma base: maintenant sur et . Maintenant c'est . Mais la méthode ne fonctionne toujours pas. vk(x)=cos2(π2h(xxi))[xi1,xi+1]i=1,,n1C1
Appliqué
La base devrait être capable de reproduire des fonctions linéaires, mais cela ne peut pas. Une fois que vous avez corrigé cela, vérifiez que les intégrales sont effectuées correctement, puis vérifiez les conditions aux limites. C1
Jed Brown
0

Clairement, le problème a un dérivé d'ordre ODD. Plus spécifiquement pour les nombres de Péclet plus importants , la matrice de rigidité peut ne pas conserver une forme `` fine '', ce qui crée des zéros lors de l'assemblage et devient donc un déterminant singulier ou parfois très petit qui est perceptible par les oscillations dans le tracé de la solution.

La solution à ce genre de problème est l'utilisation de la peine, entre autres méthodes. Plus précisément, cela s'appelle la méthode Petrov-Galerkin .

Désolé pour ma mauvaise compréhension de l'anglais.

Sohail Ahmed
la source