1D FEM per coefficiente di diffusione non lineare
Voglio risolvere con elementi finiti lineari l'equazione$$\partial_t u = \partial_{x}(a(u)\partial_xu)$$ nel dominio $t \in [0,1]$ e $x \in [-L,L]$. Qui$a(u)$ è solo una funzione di $u$.
Applicando la formulazione debole con $u(t,x)=\sum_{j} u_j(t) \varphi_j(x)$, Io ottengo $$\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$$
L'lhs non è un problema perché lo è $M \dot{U}(t)$, dove $(M)_{ij}=\int_{-L}^{L} \varphi_i(x) \varphi_j(x)dx$ e $U(t)=[u(x_1,t),\ldots,u(x_N,t)]^{T}$
Il mio grosso problema è sulle rh non so come gestire quella doppia somma in modo da avere una funzione di$U(t)$, perché ottengo un tensore $B_{ijk}=\int_{-L}^{L} \varphi_i \varphi_j \varphi_k^{'}$( c'è già stata una domanda su questo ) ma non riesco a capire come risolverlo in pratica su un computer.
Come descritto nella domanda collegata, otterrò $$M \dot{U} = (BU)U$$ma questo mi sembra solo formale, a causa di quel tensore. Qualsiasi aiuto è molto apprezzato
MODIFICA dopo la risposta knl :
@knl Ho una domanda sulla fase di ricerca della radice:
Dopo la discretizzazione del tempo, ho $u_n(x)$, quindi il problema è ancora continuo nello spazio. Dal prodotto scalare$$ (\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) $$ Voglio trovare come calcolare la soluzione tramite l'iterazione del punto fisso
Permettere $A$ la solita "matrice di rigidità" e $M$ la "matrice di massa":
$$\delta^{-1} M u_k^n + a(u_{k-1}^n) A u_k^n = \delta^{-1}M u^{n-1}$$ dove $u_k^{n}$ è il vettore dei coefficienti e $k$ è l'indice che fa riferimento all'iterazione del punto fisso.
Pertanto, trovo iterativamente $u_k^n$ risolvendo i sistemi lineari $$(\delta^{-1} M + a(u_{k-1}^n) A)u_k^n = \delta^{-1} M u^{n-1}$$
Quello che ottengo dopo l'integrazione fino al tempo $t=1$ è
che è leggermente diverso dal tuo. Non riesco a capire se c'è un errore nel mio codice, perché le iterazioni a punto fisso sembrano funzionare.
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 (ultimo) In realtà quello che ho prima è sbagliato, prima di non aggiornare la matrice nell'iterazione del punto fisso. Risolvendo questo, cioè cambiando la funzione stiffassembly, ottengo quanto segue, che è corretto, in quanto è stato controllato "tre" (con @knl fem solver, con un approccio alle differenze finite e con matematica):
Risposte
Avresti bisogno di linearizzare il problema. Preferisco farlo prima della discretizzazione ma è possibile farlo anche dopo la discretizzazione. (Sono un po 'scettico sulla linearizzazione dopo la discretizzazione perché non ho mai esaminato i dettagli. In generale, i passaggi di discretizzazione e linearizzazione non si spostano.)
Di seguito presumo che l'equazione sia effettivamente $\partial_t u = \partial_x(a(u) \partial_x u)$ e che hai la condizione al contorno $u=0$.
La forma debole è $$(\partial_t u, v) = -(a(u) \partial_x u, \partial_x v).$$Preferisco fare prima la discretizzazione del tempo in modo da vedere la struttura del problema risultante. Ad esempio, il metodo di Eulero implicito porta a$$(\delta^{-1}(u_n - u_{n-1}), v) = -(a(u_n) \partial_x u_n, \partial_x v),$$ o, equivalentemente, $$(\delta^{-1} u_n, v) + (a(u_n) \partial_x u_n, \partial_x v) = (\delta^{-1}u_{n-1}, v),$$ dove $n$ scorre nei passaggi temporali e $\delta > 0$è la dimensione del gradino. L'equazione è ancora non lineare in$u_n$e devi linearizzare. Un'opzione è eseguire un'iterazione a virgola fissa (all'interno di ogni fase temporale$n$) trovando ripetutamente $u_{k,n}$ a partire dal $$(\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),$$ dove $k$ percorre le fasi di linearizzazione e $u_{k-1,n}$è la funzione dell'iterazione precedente. Nota come ora hai due iterazioni: una per la discretizzazione temporale e una per la linearizzazione.
Ho realizzato un caso di esempio con $u(x) = \sin(\pi x)$e l'ho risolto usando il codice che conosco meglio (cioè il mio, puoi installarlo in Python usando pip install scikit-fem==2.0.0se vuoi eseguirlo):
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()
Ciò fornisce le seguenti due immagini (condizione iniziale e il risultato in $t=1$):
Ovviamente ci sono molti modi alternativi per farlo, ma questo dovrebbe darti un'idea generale.
L'equazione che hai è $$ M \dot U = B(U) U $$ dove $B$ è una matrice che dipende dalla soluzione ed è data da $$ B_{ij} = \int a(u_h) \phi_i' \phi_j' dx $$