MEF: singularidade da matriz de rigidez

{(σ^{2} (x) u^{″} (x))}^{″} = f (x), 0 ⩽ x ⩽ 1

$\left( \sigma^{2}(x) u ''(x) \right)'' = f(x), \;\;\; 0 \leqslant x \leqslant 1$

u (0) = u (1) = 0

$u(0) = u(1) = 0$

u^{″} (0) = u^{″} (1) = 0

$u''(0) = u''(1) = 0$

σ (x) ⩾ σ_{0} > 0

$\sigma(x) \geqslant \sigma_{0} > 0$

A u = f

$Au = f$

A

$A$

Seguindo o esquema do FEM, reduzo meu problema a um problema de otimização

J (u) = (A u, u) - 2 (f, u) \to min_{u}

$J(u) = (Au,u) - 2(f,u) \to \min_{u}$ Eu introduzo elementos finitos

h_{k} (x)

$h_{k}(x)$ como

v_{k} (x) = {\begin{aligned} 1 - {(\frac{x - x_{k}}{h})}^{2}, & x \in [x_{k - 1}, x_{k + 1}] \\ 0, & otherwise \end{aligned}

$v_{k}(x) = \left\{ \begin{array}{rl} 1 - \left( \frac{x-x_{k}}{h} \right)^2, & x \in [x_{k-1},x_{k+1}] \\ 0, & \text{otherwise} \end{array} \right.$ para qualquer

k = 1, \dots, n - 1

$k = 1,\ldots,n-1$ , em que

x_{k} = h k

$x_{k} = hk$ ,

h = \frac{1}{n}

$h = \frac{1}{n}$ . Os elementos finitos

v_{0} (x)

$v_{0}(x)$ e

v_{n} (x)

$v_{n}(x)$ são introduzidos de maneira semelhante.

Eu tento encontrar numericamente o vetor $\alpha$ forma que $u(x) = \sum_{k=0}^{n} \alpha_{k} v_{k}(x)$ resolva o problema de otimização. Temos

J (u) = \sum_{i = 0}^{n} \sum_{j = 0}^{n} α_{i} α_{j} (A v_{i}, v_{j}) - \sum_{i = 0}^{n} 2 α_{i} (v_{i}, f) = α^{T} V α - 2 α^{T} b \to min_{α},

$J(u) = \sum\limits_{i=0}^{n} \sum\limits_{j=0}^{n} \alpha_{i} \alpha_{j} (Av_{i},v_{j}) - \sum\limits_{i=0}^{n} 2\alpha_{i} (v_{i},f) = \alpha^{T} V \alpha - 2\alpha^{T} b \to \min\limits_{\alpha},$ onde

b_{i} = (f, v_{i})

$b_{i} = (f,v_{i})$ e

V_{i, j} = (A v_{i}, v_{j})

$V_{i,j} = (Av_{i},v_{j})$ . Após diferenciação em relação a

α

$\alpha$ , recebo

V α = b,

$V\alpha = b,$ mas aqui a matriz de rigidez

V

$V$ é singular. Então o que eu tenho que fazer? Talvez eu tenha que escolher outros elementos finitos?

finite-element Apliques
fonte

Olá, Nimza, você tem um problema de teste que sabe a solução exata? Se sim, tente resolver

V^{T} V α = V^{T} b

$V^T V \alpha = V^T b$ primeiro para testar se sua base está correta dentro do domínio, se tudo parece correto, então talvez seja o BC incorretamente posicionado que torna a matriz singular. Mas o BC parece bom para mim.

Shuhao Cao 12/04

Respostas:

Em ordem decrescente de probabilidade

Base incorreta. Na sua descrição, parece que você tem exatamente duas funções quadráticas com suporte em cada elemento. Esse espaço não é uma partição da unidade e não é (primeiras derivadas contínuas). Para discretizar seu problema de quarta ordem diretamente (em vez de reduzi-lo a um sistema de equações de segunda ordem, por exemplo), você precisará de uma base . Observe que a base deve ser capaz de reproduzir exatamente todas as funções lineares. $C^1$ $C^1$ $C^1$
Condições de contorno insuficientes. Isso será óbvio se você calcular e plotar o espaço nulo.
Montagem incorreta. Verifique o mapa dos elementos para a ordem montada para confirmar se é o que você esperava, por exemplo, que não está invertendo a orientação dos elementos.
Montagem local incorreta. Em 1D, você pode calcular analiticamente como é a matriz de rigidez do elemento (talvez para um caso simplificado) e verificar se o código a reproduz.

Jed Brown
fonte

Obrigado. 1. Acho que vou precisar de uma base porque . Então, se eu considerar apenas funções que satisfazem as condições de contorno, .

C^{2}

$C^2$

(A u, v) = \int_{0}^{1} σ^{2} (x) u^{″} (x) v^{″} (x) d x

$(Au,v) = \int_{0}^{1} \sigma^{2}(x) u''(x) v''(x) dx$

\ker A = {0}

$\ker A = \{ 0 \}$

Appliqué

Uma base é suficiente, o integrando não precisa ser contínuo. Observe que as condições de contorno nas segundas derivadas se tornarão uma integral de contorno. Você pode usar uma base para discretização direta de um problema de quarta ordem, mas precisará integrar termos de salto como nos métodos descontínuos de Galerkin para sistemas de primeira e segunda ordem. Não é um método ruim, mas é desnecessariamente complicado em 1D porque é muito fácil construir bases com qualquer ordem de continuidade (por exemplo, splines). Este artigo é um exemplo de " DG".

C^{1}

$C^1$

C^{0}

$C^0$

C^{0}

$C^0$

precisa

Está bem. minha base: agora em e . Agora é . Mas o método ainda não funciona.

v_{k} (x) = \cos^{2} (\frac{π}{2 h} (x - x_{i}))

$v_{k}(x) = \cos^2 \left(\frac{\pi}{2h}(x-x_{i}) \right)$

[x_{i - 1}, x_{i + 1}]

$[x_{i-1},x_{i+1}]$

i = 1, \dots, n - 1

$i = 1,\ldots,n-1$

C^{1}

$C^1$

Appliqué

A base deve ser capaz de reproduzir funções lineares, mas isso não pode. Depois de corrigir isso, verifique se as integrais estão sendo executadas corretamente e verifique as condições de contorno.

C^{1}

$C^1$

22812 Jed Brown

Claramente, o problema tem um derivado de ordem ímpar. Mais especificamente, para números maiores de Péclet , a matriz de rigidez pode não manter a forma 'fina', o que cria zeros durante a montagem e, portanto, fica determinante singular ou, às vezes, muito pequeno, perceptível pelas oscilações no gráfico da solução.

A solução para esse tipo de problema é o uso de penalidade, entre outros métodos. Mais especificamente, isso é chamado de método Petrov-Galerkin .

Desculpe pela minha má compreensão do inglês.

Sohail Ahmed
fonte