1D FEM pour coefficient de diffusion non linéaire
Je veux résoudre avec des éléments finis linéaires l'équation$$\partial_t u = \partial_{x}(a(u)\partial_xu)$$ dans le domaine $t \in [0,1]$ et $x \in [-L,L]$. Ici$a(u)$ est juste une fonction de $u$.
Application de la formulation faible avec $u(t,x)=\sum_{j} u_j(t) \varphi_j(x)$, J'obtiens $$\partial_t u_j(t) \int_{-L}^{L}\varphi_i(x)\varphi_j(x)dx = - \int_{-L}^{L} a\Bigl( \sum_j u_j(t) \varphi_j(x) \Bigr) \Bigl( \sum_k u_k(t) \varphi_k^{'}(x) \Bigr) \varphi_{i}^{'}(x)dx$$
Le lhs n'est pas un problème car il est $M \dot{U}(t)$, où $(M)_{ij}=\int_{-L}^{L} \varphi_i(x) \varphi_j(x)dx$ et $U(t)=[u(x_1,t),\ldots,u(x_N,t)]^{T}$
Mon gros problème est sur les rhs je ne sais pas comment gérer cette double sommation pour avoir une fonction de$U(t)$, parce que j'obtiens un tenseur $B_{ijk}=\int_{-L}^{L} \varphi_i \varphi_j \varphi_k^{'}$( il y a déjà eu une question à ce sujet ) mais je ne comprends pas comment résoudre cela en pratique sur un ordinateur.
Comme décrit dans la question liée, j'obtiendrai $$M \dot{U} = (BU)U$$mais cela me semble juste formel, à cause de ce tenseur. Toute aide est grandement appréciée
EDIT après la réponse knl :
@knl J'ai une question sur l'étape de recherche de racine:
Après la discrétisation du temps, j'ai $u_n(x)$, donc le problème est toujours continu dans l'espace. Du produit scalaire$$ (\delta^{-1} u_{k,n}, v) + (a(u_{k-1,n}) \partial_x u_{k,n}, \partial_x v) = (\delta^{-1}u_{n-1}, v) $$ Je veux trouver comment calculer la solution par itération de point fixe
Laisser $A$ la "matrice de rigidité" habituelle et $M$ la "matrice de masse":
$$\delta^{-1} M u_k^n + a(u_{k-1}^n) A u_k^n = \delta^{-1}M u^{n-1}$$ où $u_k^{n}$ est le vecteur des coefficients et $k$ est l'index faisant référence à l'itération du point fixe.
Par conséquent, je trouve itérativement $u_k^n$ en résolvant les systèmes linéaires $$(\delta^{-1} M + a(u_{k-1}^n) A)u_k^n = \delta^{-1} M u^{n-1}$$
Ce que j'obtiens après intégration jusqu'à temps $t=1$ est
ce qui est légèrement différent du vôtre. Je ne peux pas comprendre s'il y a une erreur dans mon code, car les itérations en virgule fixe semblent fonctionner.
import numpy as np
import matplotlib.pyplot as plt
def stiffassembly(M):
x = np.linspace(0,1,M+1)
diag = np.zeros(M-1) #x_1,...,x_M-1 (M-1)
subd = np.zeros(M-2)
supr = np.zeros(M-2)
h = np.diff(x)
for i in range(1,M):
diag[i-1] = 1/h[i-1] +1/h[i]
for k in range(1,M-1):
supr[k-1] = -1/h[k]
subd[k-1] = -1/h[k]
A = np.diag(subd,-1) + np.diag(diag,0) + np.diag(supr,+1)
return A
def massmatrix(N):
x = np.linspace(0,1,N+1)
diag = np.zeros(N-1) #x_1,...,x_M-1 (M-1)
subd = np.zeros(N-2)
supr = np.zeros(N-2)
h = np.diff(x)
for i in range(1,N):
diag[i-1] = (h[i-1] + h[i])/3
for k in range(1,N-1):
supr[k-1] = h[k]/6
subd[k-1] = h[k-1]/6
M = np.diag(subd,-1) + np.diag(diag,0) + np.diag(supr,+1)
return M
def inidata(x):
return np.sin(np.pi*x)
a = lambda w: (1. * w) ** 2
M = 50
x = np.linspace(0,1,M+1)
delta = 0.001
odx = 1.0/delta
tol = 1e-14
uprev = inidata(x[1:-1])
ts = 1000 #integration up to t=1.0
for n in range(ts):
print('iteration',str(n))
u = uprev.copy()
uold = u.copy() + 1
it = 0
while (np.linalg.norm(u-uold)>tol):
uold=u.copy()
u = np.linalg.solve(odx*massmatrix(M) + np.diag(a(u))@stiffassembly(M), odx*massmatrix(M)@uprev)
errnrm = np.linalg.norm(u-uold)
print(errnrm)
uprev = u.copy()
plt.plot(x,np.r_[0,u,0],'g-o',)
EDIT (dernier) En fait, ce que j'ai avant est faux, avant de ne pas mettre à jour la matrice dans l'itération du point fixe. En corrigeant cela, c'est-à-dire en changeant la fonction stiffassembly, j'obtiens ce qui suit, ce qui est juste, car il a été vérifié "trois" (avec le solveur @knl fem, avec une approche aux différences finies, et avec Mathematica):
Réponses
Vous auriez besoin de linéariser le problème. Je préfère le faire avant la discrétisation mais c'est possible de le faire aussi après la discrétisation. (Je suis un peu sceptique quant à la linéarisation après discrétisation car je n'ai jamais regardé dans les détails. En général, les étapes de discrétisation et de linéarisation ne commutent pas.)
Dans ce qui suit, je suppose que l'équation est en fait $\partial_t u = \partial_x(a(u) \partial_x u)$ et que vous avez la condition aux limites $u=0$.
La forme faible est $$(\partial_t u, v) = -(a(u) \partial_x u, \partial_x v).$$Je préfère d'abord faire la discrétisation temporelle pour que vous voyiez la structure du problème résultant. Par exemple, la méthode Euler implicite conduit à$$(\delta^{-1}(u_n - u_{n-1}), v) = -(a(u_n) \partial_x u_n, \partial_x v),$$ ou équivalent, $$(\delta^{-1} u_n, v) + (a(u_n) \partial_x u_n, \partial_x v) = (\delta^{-1}u_{n-1}, v),$$ où $n$ s'exécute sur les pas de temps et $\delta > 0$est la taille de l'étape. L'équation est toujours non linéaire dans$u_n$et vous devez linéariser. Une option consiste à faire une itération en virgule fixe (à l'intérieur de chaque pas de temps$n$) en trouvant à plusieurs reprises $u_{k,n}$ de $$(\delta^{-1} u_{k,n}, v) + (a(u_{k-1,n}) \partial_x u_{k,n}, \partial_x v) = (\delta^{-1}u_{n-1}, v),$$ où $k$ parcourt les étapes de linéarisation et $u_{k-1,n}$est la fonction de l'itération précédente. Remarquez comment vous avez maintenant deux itérations: une pour la discrétisation temporelle et une pour la linéarisation.
J'ai fait un exemple de cas avec $u(x) = \sin(\pi x)$et l'a résolu en utilisant le code que je connais le mieux (c'est-à-dire le mien, vous pouvez l'installer en Python en utilisant pip install scikit-fem==2.0.0si vous voulez l'exécuter):
from skfem import *
from skfem.helpers import *
from skfem.visuals.matplotlib import *
import numpy as np
m = MeshLine(); m.refine(5)
basis = InteriorBasis(m, ElementLineP2())
a = lambda w: (1. * w) ** 2
bilinf_stiffness = BilinearForm(lambda u, v, w: a(w['u_prev']) * dot(grad(u), grad(v)))
delta = 0.01
M = BilinearForm(lambda u, v, w: 1. / delta * u * v).assemble(basis)
load = LinearForm(lambda v, w: 1. / delta * w['u_prev'] * v)
u = project(lambda x: np.sin(np.pi * x[0]), basis_to=basis)
plot(basis, u)
for n in range(100): # 100 time steps
b = load.assemble(basis, u_prev=basis.interpolate(u))
for k in range(250): # 250 linearization loops
A = bilinf_stiffness.assemble(basis, u_prev=basis.interpolate(u))
u = solve(*condense(A + M, b, D=m.boundary_nodes()))
print("iteration {}".format(n))
plot(basis, u)
show()
Cela donne les deux images suivantes (condition initiale et résultat à $t=1$):
Il existe évidemment de nombreuses autres façons de procéder, mais cela devrait vous donner une idée générale.
L'équation que vous avez est $$ M \dot U = B(U) U $$ où $B$ est une matrice qui dépend de la solution et est donnée par $$ B_{ij} = \int a(u_h) \phi_i' \phi_j' dx $$