Finite Element Method i én dimension

Finite Element Method i én dimension#

Note

Dette projekt er lavet af Jonas Amtoft.
Med inspiration fra et FEM projekt af Martin Bendsøe. Med feedback af Jesper Bækholm Rom Poulsen.

Indledning#

Når man skal designe en bro, en bil eller en vindmøllevinge, er det afgørende at kunne forudsige, hvordan konstruktionen deformerer under belastning. Matematisk beskrives denne deformation typisk ved en differentialligning. For simple problemer kan en sådan ligning i nogle tilfælde løses analytisk, men for realistiske konstruktioner bliver problemet hurtigt for komplekst. Her bliver den endelige elementmetode, Finite Element Method (FEM), et centralt værktøj. Metoden bygger på en opdeling af domænet i et antal små elementer, hvor løsningen approksimeres lokalt. Den samlede løsning kan dermed beskrives som en sammensætning af simple funktioner, der hver er knyttet til en lille del af området. På den måde omsættes differentialligningen til et lineært ligningssystem, som kan løses numerisk. I dette projekt introduceres FEM i én dimension, hvor metoden udvikles trin for trin og implementeres i Python.

Forberedelse#

Projektet bygger på følgende emner fra Mat 1a- og Mat 1b-noterne

Mat 1a:

  • Kap. 6 - Lineære ligningssystemer og Gauss-elimination

  • Kap. 7 - Vektorer og matricer

  • Kap. 10 - Vektorrum og basis

Mat 1b:

  • Kap. 2 - Indre produktrum og ortogonalitet

  • Kap. 6 - Integration og integrationsmetoder

Projektmål#

Målet med projektet er at udvikle den endelige elementmetode (FEM) til numerisk løsning af randværdiproblemer i én dimension og undersøge metodens egenskaber gennem analyse og udvidelser.

Note

I skal skrive en sammenhængende rapport uden henvisning til opgavenumrene nedenfor.

Alle opgaver i hovedafsnittet (Opgave 1–28) er obligatoriske. I afsnittet Udvidelser vælger I selv én eller flere udvidelser at arbejde videre med, og hver udvidelse er uafhængig af de øvrige.

Lad den endelige rapport afspejle jeres interesser og ambitioner.

Randværdiproblemet#

Forestil dig en metalstang af længde 1, hvis to endepunkter holdes på samme temperatur. Langs stangen er der en fordelt varmekilde, f.eks. en elektrisk strøm der varmer materialet op. Vi ønsker at bestemme temperaturfordelingen \(u(x)\) i stangen, hvor \(x \in [0,1]\) angiver positionen. Stangens evne til at lede varme beskrives af en funktion \(a(x)\), kaldet varmeledningsevnen, og den lokale varmekilde beskrives af \(p(x)\). Under stationære betingelser er temperaturfordelingen styret af en differentialligning, der kobler \(u(x)\), \(a(x)\) og \(p(x)\).

../_images/heated-rod-schematic.png

Fig. 1 En varmeledende stang#

Et lignende eksempel er en tynd elastisk snor, fastgjort i begge ender og udsat for en fordelt belastning \(p(x)\). Her beskriver \(u(x)\) udbøjningen fra hvilestillingen, og \(a(x)\) angiver materialets stivhed. Den matematiske struktur er den samme som for varmeledningsproblemet.

../_images/elastic-string-schematic.png

Fig. 2 Nedbøjning af en tynd, elastisk snor#

Begge fysiske situationer beskrives af det samme randværdiproblem

\[ -\frac{d}{dx}\left[a(x)\,\frac{du}{dx}\right] = p(x), \qquad x \in (0,1), \quad u(0) = u(1) = 0. \tag{1} \]

Her er \(a(x)\) en differentiabel funktion med \(a(x) > 0\) for alle \(x \in [0,1]\), og \(p(x)\) er en kontinuert funktion. Randbetingelserne \(u(0) = u(1) = 0\) kaldes homogene randbetingelser.

Opgave 1#

Bestem den eksakte løsning \(u(x)\) til randværdiproblemet i ligning (1), når \(a(x)=1\) og \(p(x)=1\). Verificer, at din løsning opfylder differentialligningen og begge randbetingelser.

Opgave 2#

Skitser den eksakte løsning fra Opgave 1 i hånden for \(x \in [0,1]\). Hvor er maksimum, og hvad er \(u(1/2)\)?

Skriv derefter et kort Python-program, der beregner og plotter løsningen med np.linspace og plt.plot.

# Plot den eksakte loesning u(x) = x(1-x)/2

x = ...  # definér x-værdier med np.linspace
u = ...  # beregn den eksakte løsning
plt.figure(figsize=(6, 4))
# plt.plot, xlabel, ylabel, titel
plt.show()

Opgave 3#

Forklar kort, hvordan løsningen tolkes i henholdsvis varmeledningsmodellen og snormodellen. Hvad fortæller maksimumspunktet om den fysiske situation? Diskuter, hvad der sker, hvis kildeleddet \(p(x)\) er stærkt lokaliseret tæt ved \(x = 1\).

Opgave 4#

Betragt randværdiproblemet med \(a(x) = 1 + x\) og \(p(x) = 1\)

\[\begin{equation*} -\frac{d}{dx}\!\left[(1+x)\,\frac{du}{dx}\right] = 1, \qquad x \in (0,1), \quad u(0) = u(1) = 0. \end{equation*}\]

Find den eksakte løsning. Hvad gør det sværere end for \(a(x) = 1\)?

Den svage form#

Vi har set, at eksakte løsninger kan findes for simple \(a(x)\) og \(p(x)\). Men for komplicerede funktioner er analytiske metoder enten besværlige eller umulige. Nøglen er at omformulere randværdiproblemet. I stedet for at kræve, at differentialligningen er opfyldt punkt for punkt (den stærke form), søger vi en formulering, der kun kræver første-ordens afledede af \(u\). Strategien er at gange med en testfunktion \(v\) og integrere over \([0,1]\), derefter anvende partiel integration til at flytte én differentiering fra \(u\) over på \(v\).

Vi indfører funktionsrummet

\[ V=\{\,v:[0,1]\to\mathbb{R}\mid v \text{ er kontinuert, stykkevis } C^1,\ \text{og } v(0)=v(1)=0\,\}. \]

Det vil sige, at et element \(v\in V\) er en funktion, som er kontinuert på \([0,1]\), differentiabel på hvert delinterval, og som opfylder de homogene randbetingelser i endepunkterne.

Opgave 5#

Tag udgangspunkt i den stærke form i ligning (1). Gang begge sider med en vilkårlig testfunktion \(v \in V\) og integrer over \([0,1]\). Anvend partiel integration (se Mat 1b, Sætning 6.2.3) på venstresiden.

Konkluder, at den svage form er

\[\begin{equation*} \int_0^1 a(x)\,u'(x)\,v'(x)\,dx = \int_0^1 p(x)\,v(x)\,dx \qquad \text{for alle } v \in V. \end{equation*}\]

Vi har dermed afledt den svage form af randværdiproblemet

\[ \int_0^1 a(x)\,u'(x)\,v'(x)\,dx = \int_0^1 p(x)\,v(x)\,dx \qquad \text{for alle } v \in V. \tag{2} \]

Den svage formulering giver også mening for funktioner \(u(x)\) og \(v(x)\), som er kontinuerte og stykkevis differentiable med stykkevis kontinuerte afledede. Tilsvarende kan funktionerne \(a(x)\) og \(p(x)\) være stykkevis kontinuerte, uden at være differentiable eller overalt kontinuerte. Dermed kan formuleringen også anvendes på problemer, hvor materialet ændrer sig fra ét delinterval til et andet, fx i en sammensat stang med forskellige materialer eller tværsnit.

Opgave 6#

For \(a(x) = 1\) og \(p(x) = 1\) er den eksakte løsning \(u(x) = \frac{x(1-x)}{2}\), som vi så i Opgave 1.

(a) Verificer den svage form i ligning (2) med testfunktionen \(v(x) = x(1-x)\).

(b) Verificer med en anden testfunktion efter eget valg, f.eks. \(v(x) = \sin(\pi x)\).

(c) Argumenter for, hvorfor den svage form gælder for alle \(v \in V\), og ikke blot for de specifikke testfunktioner I har tjekket.

Opgave 7#

../_images/discontinuous-rod-schematic.png

Fig. 3 En varmeledende stang sammensat af to forskellige materialer med spring i \(a(x)\).#

Betragt en stang bestående af to materialer samlet ved \(x = \tfrac{1}{2}\) med

\[\begin{equation*} a(x) = \begin{cases} a_1, & x \in \bigl[0, \tfrac{1}{2}\bigr), \\ a_2, & x \in \bigl[\tfrac{1}{2}, 1\bigr], \end{cases} \end{equation*}\]

hvor \(a_1 \neq a_2\).

(a) Forklar, hvorfor den stærke form ikke kan fortolkes klassisk i punktet \(x = \tfrac{1}{2}\).

(b) Argumenter for, at den svage form ligning (2) er veldefineret, selv når \(a(x)\) har et spring.

Bilineær form og indre produkt#

Fra Mat 1b kender vi begrebet indre produkt (se Mat 1b, Definition 2.1.2). Vi indfører den bilineære form, dvs. en funktion af to variable, der er lineær i hver af dem separat,

\[ B(u,v) = \int_0^1 a(x)\,u'(x)\,v'(x)\,dx \tag{3} \]

samt det sædvanlige \(L^2\)-indre produkt (se Mat 1b, Eksempel 2.1.6) for reelle funktioner

\[ (f,g) = \int_0^1 f(x)\,g(x)\,dx. \]

Den svage form i ligning (2) kan nu skrives kompakt som

\[ B(u,v) = (p,v) \text{ for alle } v \in V, \quad u \in V. \tag{4} \]

Opgave 8#

(a) Vis, at \(B(\cdot,\cdot)\) defineret i ligning (3) er symmetrisk og bilineær.

(b) Vis, at \(B(v,v) \geq 0\) for alle \(v \in V\), og at \(B(v,v) = 0\) medfører \(v = 0\) i \(V\) (dvs. \(v(x) = 0\) for alle \(x\)). Konkludér at \(B(\cdot,\cdot)\) er et indre produkt på \(V\).

Hint

Til (b): \(B(v,v) = \int_0^1 a(x)[v'(x)]^2\,dx\). Hvad kan du konkludere om \(v'\), når dette integral er nul? Brug at \(v(0) = 0\).

Ligning (4) er stadig formuleret som et problem, hvor den ukendte funktion \(u\) skal findes i det uendeligt-dimensionale rum \(V\). Næste skridt er at erstatte \(V\) med et endeligt-dimensionalt underrum (se Mat 1a, Definition 10.4.1) \(V_h \subset V\), så den svage form reduceres til et lineært ligningssystem. Den metode, vi derved opnår, er Galerkin-metoden, og dens implementering med stykkevis lineære basisfunktioner er netop den endelige elementmetode (FEM).

Galerkin-metoden#

I stedet for at søge løsningen i hele \(V\) begrænser vi os til et endeligt-dimensionalt underrum \(V_h \subset V\) af dimension \(N\) med basisfunktioner \(\psi_1(x), \psi_2(x), \ldots, \psi_N(x)\) (se Mat 1a, Definition 10.2.4). Enhver funktion i \(V_h\) skrives som en linearkombination af disse.

Vi søger en approksimativ løsning på formen

\[ \hat{u}(x) = \sum_{k=1}^{N} u_k \,\psi_k(x), \qquad \hat{v}(x) = \sum_{l=1}^{N} v_l \,\psi_l(x), \tag{5} \]

hvor \(u_1, \ldots, u_N\) er ukendte koefficienter, og \(v_1, \ldots, v_N\) er vilkårlige koefficienter. Vi samler koefficienterne i vektoren \(\mathbf{u} = (u_1, \ldots, u_N)^T \in \mathbb{R}^N\).

Galerkin-betingelsen #

Vi kræver, at den svage form gælder for alle \(\hat{v} \in V_h\), dvs. for alle valg af \(v_1, \ldots, v_N\). For at \(\hat{u}\) og \(\hat{v}\) opfylder randbetingelserne, stiller vi kravet \(\psi_k(0) = \psi_k(1) = 0\) for alle \(k\).

Opgave 9#

Vis, at hvis \(\psi_k(0) = \psi_k(1) = 0\) for alle \(k = 1, \ldots, N\), så gælder \(\hat{u}(0) = \hat{u}(1) = 0\) og \(\hat{v}(0) = \hat{v}(1) = 0\) uanset valget af koefficienter.

Opgave 10#

Indsæt approksimationen i ligning (5) i den svage form, og vis, at kravet “for alle \(\hat{v} \in V_h\)” er ækvivalent med det lineære system

\[\begin{equation*} \mathbf{K}^T\mathbf{u} = \mathbf{p}, \end{equation*}\]

hvor stivhedsmatrixen \(\mathbf{K} \in \mathbb{R}^{N \times N}\) og lastvektoren \(\mathbf{p} \in \mathbb{R}^N\) er givet ved

\[\begin{equation*} K_{kl} = B(\psi_k, \psi_l) = \int_0^1 a(x)\,\psi_k'(x)\,\psi_l'(x)\,dx, \qquad k,l = 1,\ldots,N, \tag{6} \end{equation*}\]

og

\[\begin{equation*} p_l = (p, \psi_l) = \int_0^1 p(x)\,\psi_l(x)\,dx, \qquad l = 1,\ldots,N. \tag{7} \end{equation*}\]

Opgave 11#

(a) Vis, at stivhedsmatrixen \(\mathbf{K}\) er symmetrisk. Konkludér, at systemet fra Opgave 10 dermed også kan skrives som

\[\begin{equation*} \mathbf{K}\mathbf{u} = \mathbf{p}. \tag{8} \end{equation*}\]

(b) Vis, at \(\mathbf{u}^T \mathbf{K} \mathbf{u} = B(\hat{u}, \hat{u})\) for alle \(\mathbf{u} \in \mathbb{R}^N\).

(c) Vis, at \(\mathbf{K}\) er positiv definit, når \(a(x) > 0\), for enhver basis af \(V_h\). Konkludér heraf, at systemet \(\mathbf{K}\mathbf{u} = \mathbf{p}\) har en entydig løsning.

Opgave 12#

Lad \(u \in V\) være den eksakte løsning til den svage form i ligning (4), og lad \(\hat{u} \in V_h\) være Galerkin-approksimationen, dvs. løsningen \(\hat{u} \in V_h\) til ligning (4) med \(V\) erstattet af \(V_h\). Vis, at fejlen \(u - \hat{u}\) er \(B\)-ortogonal til \(V_h\), dvs.

\[\begin{equation*} B(u - \hat{u},\, v) = 0 \qquad \text{for alle } v \in V_h. \end{equation*}\]

Forklar, hvordan dette resultat kan fortolkes som en projektion i det indre produktrum \((V, B(\cdot,\cdot))\).

Et konkret eksempel med sinusbasis#

For at illustrere Galerkin-metoden vælger vi basisfunktionerne \(\psi_k(x) = \sin(k\pi x)\) for \(k = 1, \ldots, N\). Disse opfylder randbetingelserne \(\psi_k(0)=\psi_k(1)=0\) og ligger derfor i \(V\).

For \(a(x)=1\) bliver stivhedsmatricens indgange

\[ K_{kl} = \int_0^1 \psi_k'(x)\,\psi_l'(x)\,dx = \int_0^1 k\pi\cos(k\pi x)\,l\pi\cos(l\pi x)\,dx. \]

Opgave 13#

Vis, at

\[\begin{equation*} \int_0^1 \cos(k\pi x)\,\cos(l\pi x)\,dx = \begin{cases} \tfrac{1}{2} & \text{hvis } k = l, \\ 0 & \text{hvis } k \neq l. \end{cases} \end{equation*}\]

Vis derefter, at stivhedsmatrixen \(\mathbf{K}\) for sinusbasis med \(a(x)=1\) er diagonal med \(K_{kk}=\frac{k^2\pi^2}{2}\).

Opgave 14#

Cosinusfunktionerne \(\cos(k\pi x)\), \(k = 1, 2, \ldots\), er også ortogonale på \([0,1]\). Forklar, hvorfor vi ikke kan bruge \(\psi_k(x) = \cos(k\pi x)\) som basisfunktioner i Galerkin-metoden for vores randværdiproblem.

Opgave 15#

Lad \(a(x) = 1\) og \(p(x) = 1\). Beregn lastvektoren \(\mathbf{p}\) for sinusbasis med \(N = 3\), løs systemet \(\mathbf{K}\mathbf{u} = \mathbf{p}\), og sammenlign den approksimative løsning \(\hat{u}(x) = \sum_{k=1}^3 u_k \sin(k\pi x)\) med den eksakte løsning \(u(x) = \frac{x(1-x)}{2}\) i et plot.

Plottet bør vise, at Galerkin-approksimationen med blot \(N = 3\) sinusfunktioner allerede er en god tilnærmelse til den eksakte løsning. Den diagonale stivhedsmatrix gør beregningen triviel, men dette skyldes den særlige ortogonalitet af sinusfunktionerne, som gælder for \(a(x) = 1\). Hvad sker der, når \(a(x)\) varierer?

Opgave 16#

Betragt nu randværdiproblemet med stykkevis konstante funktioner

\[\begin{equation*} a(x) = \begin{cases} 2, & 0 \le x < \tfrac{1}{2}, \\ 1, & \tfrac{1}{2} \le x \le 1, \end{cases} \qquad p(x) = \begin{cases} 20, & \tfrac{2}{5} \le x < \tfrac{3}{5}, \\ 0, & \text{ellers}. \end{cases} \end{equation*}\]

(a) Beregn stivhedsmatrixen \(K_{kl} = \int_0^1 a(x)\,\psi_k'(x)\,\psi_l'(x)\,dx\) og lastvektoren \(p_k = \int_0^1 p(x)\,\sin(k\pi x)\,dx\) numerisk for sinusbasis \(\psi_k(x) = \sin(k\pi x)\) med \(N = 10\). Brug np.trapezoid (se trapezmetoden fra Opgave 5 Storedag i uge 7) til at beregne integralerne. Er \(\mathbf{K}\) stadig diagonal?

(b) Løs systemet \(\mathbf{K}\mathbf{u} = \mathbf{p}\) og plot den approksimative løsning \(\hat{u}(x) = \sum_{k=1}^N u_k \sin(k\pi x)\).

(c) Mål beregningstiden for at samle \(\mathbf{K}\) og \(\mathbf{p}\) med sinusbasis for \(N = 10, 50, 100, 200\). Hvad er \(\mathbf{K}\)’s størrelse, og hvor mange ikke-nul indgange har den?

Hint

Brug np.where til at definere de stykkevis konstante funktioner \(a(x)\) og \(p(x)\), og time modulet til tidsmåling i (c).

import time

def assemble_K_sinus(N, a_func):
    """Stivhedsmatrix for sinusbasis med generel a(x)."""
    K = np.zeros((N, N))
    x = np.linspace(0, 1, 1000)
    for k in range(1, N + 1):
        dpsi_k = ...  # den afledede af sin(k*pi*x)
        for l in range(1, N + 1):
            dpsi_l = ...  # den afledede af sin(l*pi*x)
            K[k-1, l-1] = ...  # integrer a(x)*dpsi_k*dpsi_l med np.trapezoid
    return K

def assemble_p_sinus(N, p_func):
    """Lastvektor for sinusbasis med generel p(x)."""
    x = np.linspace(0, 1, 1000)
    p_vec = np.zeros(N)
    for k in range(1, N + 1):
        psi_k = ...  # beregn psi_k = sin(k*pi*x)
        p_vec[k-1] = ...  # integrer p(x)*psi_k med np.trapezoid
    return p_vec

a_func = ...  # definer a(x) med np.where
p_func = ...  # definer p(x) med np.where

# (a) Saml K og p for N = 10
N = 10
K_sin = assemble_K_sinus(N, a_func)
p_sin = assemble_p_sinus(N, p_func)

# undersøg K, f.eks. med plt.spy eller print

# (b) Løs systemet og plot løsningen
# løs K*u = p og plot u_hat(x) = sum u_k*sin(k*pi*x)

Fra sinusbasis til lokale basisfunktioner#

I Opgave 15 viste vi, at en sinusbasis giver en diagonal \(\mathbf{K}\), når \(a(x)=1\) er konstant. I Opgave 16 så vi derimod, at denne elegante struktur forsvinder, når \(a(x)\) varierer. I det tilfælde bliver stivhedsmatrixen fuld, samlingen kræver \(N^2\) integraler, og løsningen har beregningskompleksitet \(O(N^3)\).

Begge problemer løses, hvis vi vælger basisfunktioner med kompakt support, dvs. funktioner der kun er forskellige fra nul på et lille delinterval. Så overlapper de fleste basisfunktioner ikke hinanden, og de fleste integraler \(\int a\,\psi_k'\,\psi_l'\,dx\) i stivhedsmatrixen bliver nul. Dermed bliver \(\mathbf{K}\) sparsom uanset om \(a(x)\) er konstant eller ej.

Vi ønsker derfor basisfunktioner med

  1. Kompakt support: \(\psi_k(x) = 0\) uden for et lille interval.

  2. Nemme afledede: \(\psi_k'(x)\) er stykkevis konstant.

  3. Interpolationsegenskab: \(\psi_k(x_j) = \delta_{kj}\), hvor \(\delta_{kj}\) er Kronecker-deltaet, dvs. \(\delta_{kj} = 1\) for \(k = j\) og \(\delta_{kj} = 0\) for \(k \neq j\).

Formfunktioner#

Vi vil nu vælge basisfunktioner, der giver kontinuerte og stykkevis lineære approksimative løsninger. Samtidig ønsker vi, at basisfunktionerne er lokale, så stivhedsmatrixen bliver sparsom.

Vi inddeler intervallet \([0,1]\) i \(N+1\) delintervaller med knudepunkter

\[ 0 = x_0 < x_1 < x_2 < \cdots < x_N < x_{N+1} = 1. \]

Her er \(x_0=0\) og \(x_{N+1}=1\) randpunkterne, mens de \(N\) indre knudepunkter \(x_1,\ldots,x_N\) svarer til \(N\) frihedsgrader. For et uniformt net er \(x_k=kh\) med \(h=1/(N+1)\) (\(h\) er den konstante afstand mellem knudepunkterne), men definitionen gælder også for vilkårlige knudepunkter.

Til hvert indre knudepunkt \(x_k\) knytter vi en formfunktion \(\phi_k(x)\), \(k=1,\ldots,N\), som er kontinuert, stykkevis lineær og kun forskellig fra nul på intervallet \([x_{k-1},x_{k+1}]\). Den opfylder

\[ \phi_k(x_{k-1}) = 0, \qquad \phi_k(x_k) = 1, \qquad \phi_k(x_{k+1}) = 0, \]

og er lineær på intervallerne \([x_{k-1},x_k]\) og \([x_k,x_{k+1}]\).

../_images/hat-functions-schematic.png

Fig. 4 Venstre: De tre formfunktioner \(\phi_1, \phi_2, \phi_3\) for \(N = 3\). Højre: En stykkevis lineær approksimation \(\hat{u}(x) = \sum_{k=1}^3 u_k\,\phi_k(x)\).#

Opgave 17#

(a) Vis, at \(\phi_k\) er en stykkevis lineær og kontinuert funktion givet ved

\[\begin{equation*} > \phi_k(x) = \begin{cases} > \dfrac{x - x_{k-1}}{x_k - x_{k-1}} & x_{k-1} \le x < x_k, \\[6pt] > \dfrac{x_{k+1} - x}{x_{k+1} - x_k} & x_k \le x < x_{k+1}, \\[6pt] > 0 & \text{ellers.} > \end{cases} > \qquad \text{(9)} > \end{equation*}\]

(b) Vis, at hvis \(\hat{u}(x) = \sum_{k=1}^N u_k\,\phi_k(x)\), da er \(\hat{u}\) kontinuert, stykkevis lineær og opfylder \(\hat{u}(0) = 0\) og \(\hat{u}(1) = 0\). Vis desuden, at \(\hat{u}(x_k) = u_k\) for \(k = 1, \ldots, N\).

(c) Beregn den afledede \(\phi_k'(x)\) på hvert delinterval. Hvad bliver \(\phi_k'(x)\) for et uniformt net med \(h = x_k - x_{k-1}\)?

Når vi bruger \(\phi_1, \ldots, \phi_N\) som basisfunktioner i Galerkin-metoden, får vi FEM-approksimationen \(\hat{u}(x) = \sum_{k=1}^N u_k\,\phi_k(x)\), hvor koefficienterne bestemmes ved \(\mathbf{K}\mathbf{u} = \mathbf{p}\). I det følgende arbejder vi med et uniformt net, \(h = 1/(N+1)\) og \(x_k = kh\).

Opgave 18#

Implementer og plot de \(N = 3\) formfunktioner \(\phi_1, \phi_2, \phi_3\) på intervallet \([0,1]\) for det uniforme net.

N = 3
h = 1.0 / (N + 1)
x_plot = np.linspace(0, 1, 500)

plt.figure(figsize=(8, 4))
for k in range(1, N + 1):
    x_k = ...  # knudepunktets position
    phi_k = ...  # formfunktionen phi_k
    plt.plot(x_plot, phi_k, label=f'$\\phi_{k}$')

plt.xlabel('$x$')
plt.ylabel('$\\phi_k(x)$')
plt.title(f'Formfunktioner for $N = {N}$')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Opgave 19#

Vis, at \(\phi_k\) og \(\phi_l\) har overlappende kompakt support hvis og kun hvis \(|k - l| \leq 1\), og at \(K_{kl} = 0\) når \(|k - l| \geq 2\). Hvilken matrixstruktur har \(\mathbf{K}\)?

Sammensætning og løsning#

Fra Opgave 19 ved vi, at \(\mathbf{K}\) er tridiagonal. Nu beregner vi de ikke-nul indgange eksplicit for et uniformt net med konstant \(a(x) = a\). Idéen er at gå element for element. For hvert delelement \([x_k, x_{k+1}]\) er der præcis to formfunktioner med ikke-nul afledet, og deres bidrag samles i den globale matrix \(\mathbf{K}\).

Opgave 20#

Benyt \(\phi_k'(x)\) fra Opgave 17 til at beregne de tre typer af indgange i stivhedsmatrixen for konstant \(a(x) = a\):

(a) Diagonalindgangen \(K_{kk} = \int_0^1 a\,(\phi_k')^2\,dx\).

(b) Nabo-indgangen \(K_{k,k+1} = \int_0^1 a\,\phi_k'\,\phi_{k+1}'\,dx\) (og vis \(K_{k+1,k} = K_{k,k+1}\)).

(c) Skriv \(\mathbf{K}\) eksplicit op for \(N = 3\) med \(a = 1\) og \(h = 1/4\).

Opgave 21#

Beregn lastvektoren for \(p(x) = 1\):

\[\begin{equation*} p_k = \int_0^1 \phi_k(x)\,dx = \int_{x_{k-1}}^{x_{k+1}} \phi_k(x)\,dx. \end{equation*}\]

Hvad er \(\mathbf{p}\) for \(N = 3\)?

Opgave 22#

(a) Implementer funktionen assemble_K(N, a=1.0), der definerer stivhedsmatrixen \(\mathbf{K}\) for \(N\) indre knudepunkter og konstant \(a\).

(b) Implementer funktionen assemble_p(N), der definerer lastvektoren \(\mathbf{p}\) for \(p(x) = 1\).

(c) Løs FEM-systemet for \(N = 5\) med np.linalg.solve og plot FEM-løsningen mod den eksakte løsning \(u(x) = \frac{x(1-x)}{2}\).

def assemble_K(N, a=1.0):
    """Saml stivhedsmatricen K for N indre knudepunkter, konstant a."""
    h = 1.0 / (N + 1)
    K = np.zeros((N, N))
    for k in range(N):
        K[k, k] = ...  # diagonalindgang
        if k > 0:
            K[k, k-1] = ...  # nedre naboindgang
        if k < N - 1:
            K[k, k+1] = ...  # øvre naboindgang
    return K

def assemble_p(N):
    """Saml lastvektoren p for p(x) = 1."""
    h = ...  # skridtlængde
    return ...  # lastvektor med alle indgange lig h

# Løs og plot for N = 5
N = 5

K = ...  # saml stivhedsmatrixen
p_vec = ...  # saml lastvektoren
u = np.linalg.solve(K, p_vec)

# plot FEM-løsningen sammen med den eksakte løsning

Opgave 23#

Lav et plot for hver af værdierne \(N=5\), \(N=10\) og \(N=20\), og vis i hvert tilfælde både FEM-løsningen og den eksakte løsning.

# plot FEM-løsninger for N = 5, 10, 20
# sammen med den eksakte løsning u(x) = x(1-x)/2

FEM-approksimationen forbedres tydeligt med voksende \(N\). Kurverne nærmer sig den eksakte løsning, og allerede for \(N = 20\) er forskellen svær at se med det blotte øje.

Generalisering til varierende funktioner#

Vi har hidtil brugt konstant \(a(x) = a\) og kildefunktionen \(p(x) = 1\), hvor integralerne i \(\mathbf{K}\) og \(\mathbf{p}\) kunne beregnes analytisk. Men randværdiproblemet i ligning (1) gælder for generelle \(a(x)\) og \(p(x)\), og mange fysiske situationer kræver netop denne generalitet, for eksempel stangen med varierende materialeegenskaber fra Opgave 7.

For generelle funktioner \(a(x)\) og \(p(x)\) kan integralerne i \(K_{kl}\) og \(p_k\) sjældent beregnes analytisk. Vi approksimerer dem i stedet. Da \(k\) og \(l\) nu bruges som knudeindeks i \(K_{kl}\), skifter vi til \(e\) som elementindeks for at undgå forveksling. Da \(\phi_k'\) er konstant på hvert element, er det tilstrækkeligt at approksimere \(a(x)\) med dens værdi i elementets midtpunkt \(\bar{x}_e = (x_e + x_{e+1})/2\)

\[ \int_{x_e}^{x_{e+1}} a(x)\,\phi_k'(x)\,\phi_l'(x)\,dx \approx a(\bar{x}_e)\cdot \phi_k' \cdot \phi_l' \cdot h_e, \]

hvor \(h_e = x_{e+1} - x_e\) er elementbredden. Denne fremgangsmåde generaliserer samtidig til ikke-uniforme net, hvor hvert element kan have sin egen bredde \(h_e\). Det giver mulighed for at placere flere knudepunkter, hvor løsningen varierer hurtigt, og færre, hvor den er glat.

Important

Når \(a(x)\) har et spring ved et punkt \(x = c\), bør nettet have en knude i springpunktet. Ellers spænder ét element over to materialer, og midtpunktapproksimationen repræsenterer ikke funktionen korrekt på dette element.

Opgave 24#

Implementer funktionen assemble_K_general(nodes, a_func), der samler stivhedsmatrixen for et vilkårligt net nodes \(= [x_0, x_1, \ldots, x_{N+1}]\) og en funktion a_func. Brug midtpunktreglen til at evaluere \(a(x)\) på hvert element.

def assemble_K_general(nodes, a_func=None):
    """Saml K for vilkårligt net og funktion a(x)."""
    N = len(nodes) - 2
    K = np.zeros((N, N))
    if a_func is None:
        a_func = lambda x: 1.0
    for e in range(N + 1):
        h_e = ...  # elementbredde
        x_mid = ...  # midtpunkt af elementet
        a_e = ...  # a(x) evalueret i midtpunktet
        i = e - 1
        j = e
        # opdater K[i,i], K[j,j], K[i,j], K[j,i]
        # Husk at tjekke om i og j er indre knudepunkter
    return K

Opgave 25#

Implementer assemble_p_general(nodes, p_func), der samler lastvektoren med numerisk integration (brug np.trapezoid) for en generel kildefunktion \(p(x)\).

def assemble_p_general(nodes, p_func=None):
    """Saml lastvektoren p for vilkårligt net og kildefunktion p(x)."""
    N = len(nodes) - 2
    p_vec = np.zeros(N)
    if p_func is None:
        p_func = lambda x: 1.0
    for k in range(N):
        x_left = ...  # venstre endepunkt af phi_k's support
        x_right = ...  # højre endepunkt af phi_k's support
        x_int = np.linspace(x_left, x_right, 200)
        x_k = nodes[k + 1]
        h_left = ...  # bredde af venstre delelement
        h_right = ...  # bredde af højre delelement
        phi_k = np.where(
            x_int <= x_k,
            (x_int - x_left) / h_left,
            (x_right - x_int) / h_right
        )
        p_vec[k] = ...  # integrer p(x)*phi_k med np.trapezoid
    return p_vec

Vi har nu både en Galerkin-implementering med sinusbasis og en generel FEM-samling med lokale formfunktioner. Først på dette tidspunkt er det derfor rimeligt at sammenligne beregningsomkostningerne for de to metoder på det samme randværdiproblem med varierende koefficienter.

Opgave 26#

I har nu implementeret både Galerkin-metoden med sinusbasis og en generel FEM-samling med lokale formfunktioner.

Betragt randværdiproblemet med den glatte variable koefficient

\[\begin{equation*} >a(x)=1+x, \qquad p(x)=1. >\end{equation*}\]

(a) Mål beregningstiden for at samle og løse Galerkin-systemet med sinusbasis for \(N = 10, 50, 100, 200\).

(b) Mål beregningstiden for at samle og løse FEM-systemet med lokale formfunktioner for de samme værdier af \(N\), brug et uniformt net.

(c) Sammenlign størrelsen og strukturen af stivhedsmatricerne i de to tilfælde. Hvor mange ikke-nul indgange har de?

(d) Forklar forskellen i beregningstid ud fra basisfunktionernes support og den resulterende matrixstruktur.

import time

a_func = lambda x: 1 + x
p_func = lambda x: 1 + 0*x

print(f"{'N':>6}  {'Sinusbasis (s)':>16}  {'Formfkt. (s)':>14}")
print("-" * 42)
for N in [10, 50, 100, 200]:
    # Sinusbasis
    t0 = time.time()
    K_sin = ...   # saml med assemble_K_sinus
    p_sin = ...   # saml med assemble_p_sinus
    u_sin = ...   # løs systemet
    t_sin = time.time() - t0

    # Formfunktioner på uniformt net
    nodes = ...   # uniformt net med N indre knudepunkter
    t0 = time.time()
    K_fem = ...   # saml med assemble_K_general
    p_fem = ...   # saml med assemble_p_general
    u_fem = ...   # løs systemet
    t_fem = time.time() - t0

    print(f"{N:>6}  {t_sin:>16.2e}  {t_fem:>14.2e}")

Opgave 27#

Betragt den stykkevis konstante funktion

\[\begin{equation*} a(x) = \begin{cases} 2, & 0 \le x \le \tfrac{1}{2}, \\ 1, & \tfrac{1}{2} < x \le 1. \end{cases} \end{equation*}\]

Konstruér et net med knudepunktet \(x = 0.5\) og løs randværdiproblemet med \(p(x) = 1\). Plot FEM-løsningen og sammenlign med konstant \(a(x) = 1\) og \(a(x) = 2\).

Opgave 28#

Betragt kildefunktionen \(p(x) = \exp\!\bigl(-\frac{(x-0.5)^2}{2 \cdot 0.05^2}\bigr)\), en smal Gauss-klokke, med \(a(x) = 1\). Den eksakte løsning til dette randværdiproblem er implementeret i funktionen u_exact_gaussian nedenfor.

(a) Løs med et uniformt net med \(N = 8\). Plot FEM-løsningen sammen med den eksakte løsning.

(b) Konstruér et ikke-uniformt net med \(N = 8\), hvor knuderne er tættere placeret nær \(x = 0.5\) og mere spredte mod endepunkterne. Plot FEM-løsningen sammen med den eksakte løsning.

(c) Forklar, hvorfor det ikke-uniforme net giver bedre nøjagtighed med det samme antal knudepunkter.

from scipy.special import erf

def u_exact_gaussian(x, sigma=0.05):
    x = np.asarray(x)
    z = (x - 0.5) / (np.sqrt(2) * sigma)
    c = 1 / (2 * np.sqrt(2) * sigma)

    return (
        sigma**2 * (
            np.exp(-1 / (8 * sigma**2))
            - np.exp(-((x - 0.5)**2) / (2 * sigma**2))
        )
        + sigma * np.sqrt(np.pi / 2) * (
            0.5 * erf(c) - (x - 0.5) * erf(z)
        )
    )

I har nu udvidet FEM-metoden fra det oprindelige randværdiproblem med konstante funktioner til en generel løser, der kan håndtere varierende materialeegenskaber og ikke-uniforme net. Denne fremgangsmåde med svag form, Galerkin-diskretisering, elementvis samling og numerisk løsning er kernen i FEM og genfindes i to og tre dimensioner. I de følgende udvidelser kan I arbejde videre med aspekter, som metoden bygger på.

Udvidelser#

Herunder beskrives en række uafhængige udvidelser af projektet. I vælger selv, hvilke der interesserer jer, og ingen af dem kræver, at I har lavet de øvrige.

Udvidelse A - Fejlanalyse#

Når man siger, at man har fundet en “approksimation” til løsningen af et randværdiproblem, underforstår man, at approksimationen ligger “tæt” på den rigtige løsning \(u\). Men udtrykket “tæt på” skal præciseres. Vi indfører \(L^2\)-fejlnormen

\[ F(N) = \left[\int_0^1 \bigl(u(x) - \hat{u}(x)\bigr)^2\,dx\right]^{1/2}, \tag{10} \]

der måler den samlede afvigelse mellem den eksakte løsning \(u\) og FEM-approksimationen \(\hat{u}\) over hele \([0,1]\). Tallet \(F(N)\) afhænger af antallet af indre knudepunkter \(N\), og vi forventer, at \(F(N) \to 0\) for \(N \to \infty\). Spørgsmålet er hvor hurtigt.

I denne udvidelse undersøger vi \(F(N)\) først numerisk for det specifikke problem \(a(x) = 1\), \(p(x) = 1\) (hvor vi kender den eksakte løsning fra Opgave 1), og derefter udleder vi konvergensraten analytisk for dette problem.

Opgave A.1#

Sæt \(a(x) = 1\) og \(p(x) = 1\). For \(N = 2\) er \(h = 1/3\) med knudepunkter \(x_1 = 1/3\) og \(x_2 = 2/3\).

(a) Opstil \(\mathbf{K}\) og \(\mathbf{p}\) for \(N = 2\), løs systemet \(\mathbf{K}\mathbf{u} = \mathbf{p}\), og vis at FEM-løsningens nodepunktsværdier er \(\hat{u}(x_1) = \hat{u}(x_2) = 1/9\).

(b) Beregn \(F(2)\) analytisk ved at integrere element for element. Udnyt symmetrien i problemet.

Opgave A.2#

Vi fortsætter med problemet hvor \(a(x)=1\) og \(p(x)=1\).

(a) Skriv en funktion fem_solution(x_fine, u_coeffs, N), der evaluerer FEM-løsningen på et net ved stykkevis lineær interpolation med np.interp.

(b) Skriv en funktion l2_error(N), der beregner \(F(N)\) numerisk med np.trapezoid.

(c) Beregn \(F(N)\) for \(N = 4, 8, 16, 32, 64, 128, 256\). Lav et log-log-plot af \(F(N)\) som funktion af \(N\), tilpas en ret linje med np.polyfit og angiv hældningen. Vis ved numeriske eksperimenter, at \(F(N) \propto N^{-2}\).

def fem_solution(x_fine, u_coeffs, N):
    """Evaluer FEM-løsningen på et fint gitter via stykkevis lineær interpolation."""
    u_full = np.concatenate([[0.0], u_coeffs, [0.0]])
    x_nodes = np.linspace(0, 1, N + 2)
    return np.interp(x_fine, x_nodes, u_full)

def l2_error(N):
    """Beregn L2-fejlnormen F(N) = ||u - u_hat||_2."""
    K = ...  # saml stivhedsmatrixen
    p_vec = ...  # saml lastvektoren
    u = ...  # løs systemet
    x_fine = ...  # fint gitter til integration
    u_approx = ...  # evaluer FEM-løsningen på x_fine
    u_exact = ...  # eksakt løsning u(x) = x(1-x)/2
    return ...  # returner L2-normen af fejlen (brug np.trapezoid)

Hældningen \(\approx -2\) bekræfter \(F(N) \propto N^{-2}\).

Opgave A.3#

For et uniformt net viser de numeriske eksperimenter, at \(F(N)\propto N^{-2}\) for \(a(x)=1\) og \(p(x)=1\). Vis dette analytisk ved at beregne \(F(N)\) eksakt og vis, at

\[\begin{equation*} F(N)=\frac{C}{(N+1)^2} \end{equation*}\]

for en passende konstant \(C\).

Vi har nu set numerisk og analytisk, at \(F(N) \propto N^{-2}\). Men vi har ikke overvejet, om \(L^2\)-normen er det eneste rimelige fejlmål. Det gør vi i den afsluttende opgave.

Opgave A.4#

Vi vælger som nævnt at måle fejlen med \(L^2\)-normen \(F(N) = \bigl[\int_0^1 (u - \hat{u})^2\,dx\bigr]^{1/2}\). Et andet naturligt valg er \(L^\infty\)-normen (supremumsnormen)

\[\begin{equation*} G(N) = \max_{x \in [0,1]} |u(x) - \hat{u}(x)|. \end{equation*}\]

(a) Forklar med egne ord forskellen mellem de to normer. Hvad måler \(F(N)\), og hvad måler \(G(N)\)? Giv et eksempel, hvor de to normer giver vidt forskellige vurderinger af fejlen.

(b) For problemet \(-u'' = 1\) med \(u(0)=u(1)=0\) rammer FEM de eksakte nodepunktsværdier. Hvad er \(G(N)\) i dette tilfælde, og hvad er \(F(N)\)? Giver de to normer den samme konvergensrate?

(c) Argumenter for, at \(L^2\)-normen er et mere retfærdigt mål for den samlede approksimationskvalitet end \(L^\infty\)-normen. Hvornår ville \(L^\infty\) være at foretrække?

Udvidelse B - Tridiagonal struktur og effektiv løsning#

I Opgave 19 viste I, at stivhedsmatrixen \(\mathbf{K}\) er tridiagonal, dvs. kun tre af matrixens diagonaler er ikke-nul. Denne struktur er en direkte konsekvens af, at formfunktionerne \(\phi_k\) har kompakt support og kun overlapper med deres nærmeste naboer. I Opgave 26 observerede I, at np.linalg.solve er langt hurtigere for den tridiagonale \(\mathbf{K}\) end for den fulde sinusbasis-matrix. I denne udvidelse undersøger vi, hvorfor dette er tilfældet, og implementerer en specialiseret algoritme, der udnytter den tridiagonale struktur direkte.

For en fuld \(N \times N\)-matrix koster Gauss-elimination \(O(N^3)\) operationer. Men for en tridiagonal matrix kan Thomas-algoritmen (en specialiseret Gauss-elimination) løse systemet i blot \(O(N)\) operationer. Se Thomas-algoritmen for en detaljeret beskrivelse.

plt.figure(figsize=(5, 5))
plt.spy(assemble_K(10), markersize=8)
plt.title('Ikke-nul indgange i $\\mathbf{K}$ ($N=10$)')
plt.tight_layout()
plt.show()

Opgave B.1#

(a) Vis, at en \(N \times N\) tridiagonal matrix har \(3N - 2\) ikke-nul indgange. Sammenlign med \(N^2\) for en fuld matrix.

(b) For \(N = 1000\), hvad er forholdet \(N^2 / (3N-2)\)?

Thomas-algoritmen består af to trin

  1. Fremad-eliminering: Eliminér subdiagonalen.

  2. Bagud-substitution: Løs det øvre-triangulære system bagfra.

Hvert trin kræver \(O(N)\) operationer.

Opgave B.2#

Implementer Thomas-algoritmen og verificer den mod np.linalg.solve.

def thomas(d, u_diag, l_diag, b):
    """
    Løs det tridiagonale system Ax = b med Thomas-algoritmen.

    Parametre
    ---------
    d      : ndarray, form (N,)    - Hoveddiagonal
    u_diag : ndarray, form (N-1,) - Superdiagonal (øvre)
    l_diag : ndarray, form (N-1,) - Subdiagonal (nedre)
    b      : ndarray, form (N,)    - Højreside

    Returnerer
    ----------
    x : ndarray, form (N,) - Løsningsvektoren
    """
    N = len(d)
    d = d.copy()
    b = b.copy()

    # Fremad-eliminering
    for i in range(1, N):
        factor = ...  # eliminationsfaktor l_diag[i-1] / d[i-1]
        d[i] = ...  # opdater diagonal
        b[i] = ...  # opdater højreside

    # Bagud-substitution
    x = np.zeros(N)
    x[-1] = ...  # sidste komponent
    for i in range(N-2, -1, -1):
        x[i] = ...  # bagud-substitution

    return x

# Verificer din Thomas-algoritme mod np.linalg.solve

Opgave B.3#

Undersøg egenskaberne ved stivhedsmatrixen \(\mathbf{K}\) og overvej, om Thomas-algoritmen fra Opgave B.2 kan forenkles. Implementer en optimeret version og verificer den mod np.linalg.solve.

Note

I praksis kan man bruge scipy.linalg.solve_banded, der implementerer en optimeret Thomas-algoritme i Fortran (LAPACK). Se SciPy-dokumentationen.

Udvidelse C - Fjedermodellen#

FEM-stivhedsmatrixen \(\mathbf{K}\) har en karakteristisk tridiagonal form, hvor hvert indre knudepunkt kun kobler til sine to naboer. Denne struktur er ikke tilfældig, men afspejler at formfunktionerne \(\phi_k\) kun overlapper med deres nærmeste naboer. Det viser sig, at et rent diskret mekanisk system bestående af en kæde af fjedre giver præcis den samme matrixstruktur. I denne udvidelse undersøger vi denne sammenhæng.

Betragt en kæde af \(N+1\) identiske fjedre med fjederkonstant \(\gamma > 0\), spændt vandret ud mellem to faste vægge. De \(N\) indre knuder \(y_1, \ldots, y_N\) kan forskydes, mens endepunkterne er fastholdte: \(y_0 = y_{N+1} = 0\).

../_images/spring-chain-schematic.png

Fig. 5 Fjedre på række#

Ved hver indre knude virker en ydre kraft \(F_k\). Tænk på det som at nogen skubber med en finger på knude \(k\) - det er en kraft pålagt udefra, i modsætning til de indre fjederkræfter, som fjedrene selv udøver når de strækkes eller komprimeres.

Når systemet er i ligevægt, skal summen af alle kræfter på hver knude være nul. Dette princip kaldes en kraftbalance. For knude \(k\) bidrager de to tilstødende fjedre hver med en kraft ifølge Hookes lov:

  • Venstre fjeder (mellem knude \(k{-}1\) og \(k\)): trækker med kraften \(\gamma(y_k - y_{k-1})\)

  • Højre fjeder (mellem knude \(k\) og \(k{+}1\)): trækker med kraften \(\gamma(y_k - y_{k+1})\)

Den samlede fjederkraft på knude \(k\) er summen af disse to bidrag:

\[ \gamma(y_k - y_{k-1}) + \gamma(y_k - y_{k+1}) = -\gamma\,y_{k-1} + 2\gamma\,y_k - \gamma\,y_{k+1}. \]

Kraftbalancen siger, at denne samlede fjederkraft skal balancere den ydre kraft \(F_k\):

\[ -\gamma\,y_{k-1} + 2\gamma\,y_k - \gamma\,y_{k+1} = F_k, \qquad k = 1, \ldots, N. \tag{11} \]

Opgave C.1#

(a) Skriv systemet i ligning (11) på matrixform \(\mathbf{K}_{\text{fjeder}}\,\mathbf{y} = \mathbf{F}\) for \(N = 3\).

(b) Sammenlign \(\mathbf{K}_{\text{fjeder}}\) med FEM-stivhedsmatrixen \(\mathbf{K}\) fra Opgave 20 for \(a = 1\). Hvilken værdi af \(\gamma\) giver identiske matricer? Forklar, hvorfor de to matricer har samme struktur.

Opgave C.2#

Løs fjedersystemet numerisk for \(N = 7\) og \(\gamma = 10\). Vælg kraftfordelingen \(F_k = h\) (med \(h = 1/(N+1)\)), så højresiden matcher FEM-lastvektoren. Bestem den værdi af \(a\), der giver \(\mathbf{K}_{\text{fjeder}} = \mathbf{K}_{\text{FEM}}\) (brug resultatet fra Opgave C.1), og plot forskydningerne \(y_1, \ldots, y_N\) sammen med FEM-løsningen.

Udvidelse D - Bjælkeligningen#

I den obligatoriske del har vi arbejdet med andenordens randværdiproblemer, hvor den svage form kun involverer første afledede. De stykkevis lineære formfunktioner \(\phi_k\) er tilstrækkelige, fordi \(\phi_k'\) er veldefineret.

Vi vil nu undersøge fjerdeordens randværdiproblemer, som blandt andet bruges til at beskrive bjælkers nedbøjning. Her involverer den svage form anden-afledede, og basisfunktionerne skal derfor have en veldefineret og kontinuert førstafledet, dvs. de skal være \(C^1\), ikke blot \(C^0\). De stykkevis lineære formfunktioner kan ikke bruges, fordi deres afledede har spring i knudepunkterne. I stedet anvender vi Hermite-kubiske basisfunktioner, som er stykkevis kubiske (tredjegrads) polynomier.

Bjælkeproblemet#

En bjælke af længde 1, fastspændt i begge ender, udsættes for en fordelt belastning \(P(x)\). Nedbøjningen \(w(x)\) opfylder det fjerdeordens randværdiproblem

\[ \frac{d^2}{dx^2}\left[EI(x)\,\frac{d^2 w}{dx^2}\right] = P(x), \qquad x \in (0,1), \tag{12} \]

hvor \(EI(x) > 0\) er bjælkens bøjningsstivhed. Når \(EI\) er konstant, reducerer ligningen til \(EI\,w''''(x) = P(x)\). De fire randbetingelser

\[ w(0) = 0, \quad w'(0) = 0, \quad w(1) = 0, \quad w'(1) = 0 \]

udtrykker, at bjælken hverken kan forskydes eller rotere ved endepunkterne. Bemærk, at et fjerdeordens problem kræver fire randbetingelser mod to for andenordens problemer.

Opgave D.1#

Vis, at funktionen \(w(x) = \frac{x^4}{24} - \frac{x^3}{12} + \frac{x^2}{24}\) løser ligning (12) med \(EI = 1\) og \(P(x) = 1\). Verificer alle fire randbetingelser.

Opgave D.2#

Indfør funktionsrummet \(W = \{\,v : [0,1]\to\mathbb{R} \mid v \text{ er } C^1, \text{ stykkevis } C^2,\; v(0)=v'(0)=v(1)=v'(1)=0\,\}\). Udled den svage form for bjælkeligningen ved at gange med \(v \in W\) og integrere over \([0,1]\).

Hermite-kubiske basisfunktioner#

Den svage form involverer \(w''\) og \(v''\), som skal være veldefinerede og integrable. Det kræver, at basisfunktionerne har kontinuerte første afledede (\(C^1\)), og de stykkevis lineære formfunktioner \(\phi_k\) fra hoveddelen er kun \(C^0\) og kan derfor ikke bruges.

Vi indfører et uniformt net med \(N\) indre knudepunkter, skridtlængde \(h = 1/(N+1)\) og knudepunkter \(x_k = kh\) for \(k = 1, \ldots, N\). For hvert knudepunkt definerer vi to basisfunktioner \(\zeta_k\) og \(\eta_k\).

\[\begin{split} \zeta_k(x) = \begin{cases} \left(\left|\dfrac{x - x_k}{h}\right| - 1\right)^{\!2} \left(2\left|\dfrac{x - x_k}{h}\right| + 1\right), & x_{k-1} \leq x \leq x_{k+1}, \\[6pt] 0, & \text{ellers,} \end{cases} \end{split}\]
\[\begin{split} \eta_k(x) = \begin{cases} (x - x_k)\left(\left|\dfrac{x - x_k}{h}\right| - 1\right)^{\!2}, & x_{k-1} \leq x \leq x_{k+1}, \\[6pt] 0, & \text{ellers.} \end{cases} \end{split}\]

Begge funktioner er nul uden for \([x_{k-1}, x_{k+1}]\) og er tredjegradspolynomier på hvert af delintervallerne \([x_{k-1}, x_k]\) og \([x_k, x_{k+1}]\).

Opgave D.3#

(a) Plot formfunktionerne \(\zeta_k(x)\) og \(\eta_k(x)\) for ét indre knudepunkt \(x_k\). Vis, at de er differentiable funktioner med kontinuert afledet.

(b) Aflæs fra plottet og verificer analytisk, at

\[\zeta_k(x_k) = 1, \quad \zeta_k'(x_k) = 0, \quad \eta_k(x_k) = 0, \quad \eta_k'(x_k) = 1.\]

Vi følger nu samme fremgangsmåde som i hoveddelen ved at skrive den ukendte funktion som en linearkombination af basisfunktioner og indsætte dette udtryk i den svage form, hvilket igen fører til et lineært system. Forskellen er, at vi nu har to basisfunktioner per knudepunkt, \(\zeta_k\) og \(\eta_k\), i stedet for én, \(\phi_k\), og dermed dobbelt så mange ubekendte. Bemærk, at vi kun definerede \(\zeta_k\) og \(\eta_k\) for de \(N\) indre knudepunkter \(x_1, \ldots, x_N\). Ved randknuderne \(x_0\) og \(x_{N+1}\) findes der tilsvarende lokale Hermite-former, men de tilhørende frihedsgrader bestemmes entydigt af randbetingelserne \(w(0)=w'(0)=w(1)=w'(1)=0\) og indgår derfor ikke blandt de globale ubekendte.

Lad \(\mathbf{w} = (w_1, w_2, \ldots, w_{2N})\) være \(2N\) reelle koefficienter. Approksimationen skrives som

\[ \hat{w}(x) = \sum_{k=1}^{N} w_{2k-1}\,\zeta_k(x) + \sum_{k=1}^{N} w_{2k}\,\eta_k(x). \]

Indsættes approksimationen \(\hat{w}\) i den svage form fra Opgave D.2 og vælges testfunktionerne \(v\) som basisfunktionerne selv, fås et \(2N \times 2N\) lineært system

\[ \mathbf{K}\mathbf{w} = \mathbf{f}, \]

hvor \(\mathbf{f}\) er lastvektoren med indgange \(f_m = \int_0^1 P(x)\,\psi_m(x)\,dx\) og \(\psi_m\) løber over basisfunktionerne \(\zeta_1, \eta_1, \zeta_2, \eta_2, \ldots, \zeta_N, \eta_N\).

Stivhedsmatrixen \(\mathbf{K}\) har \(2N \times 2N\) indgange. Hver indgang fås ved at beregne et integral mellem to af de \(2N\) basisfunktioner. Denne opbygning er helt parallel til hoveddelen - den eneste forskel er, at integralerne nu involverer anden-afledede i stedet for første afledede.

Ligesom i hoveddelen samles \(\mathbf{K}\) ved at beregne bidrag fra hvert delelement og lægge dem sammen. På hvert delinterval \([x_k, x_{k+1}]\) er der fire lokale basisfunktioner (\(\zeta_k\), \(\eta_k\), \(\zeta_{k+1}\), \(\eta_{k+1}\)), som er de eneste basisfunktioner, der er forskellige fra nul på dette interval. Elementstivhedsmatrixen - dvs. \(4 \times 4\)-matrixen af integraler over netop dette ene delinterval - er

\[\begin{equation*} \mathbf{K}^{(k)} = \frac{EI}{h^3} \begin{pmatrix} 12 & 6h & -12 & 6h \\ 6h & 4h^2 & -6h & 2h^2 \\ -12 & -6h & 12 & -6h \\ 6h & 2h^2 & -6h & 4h^2 \end{pmatrix} \qquad \text{(13)} \end{equation*}\]

med frihedsgradsordening \([w_{2k-1},\; w_{2k},\; w_{2k+1},\; w_{2k+2}]\), hvor \(EI\) er antaget konstant på elementet.

Opgave D.4#

Indgang \((1,1)\) i elementstivhedsmatrixen (13) svarer til basisfunktionen \(\zeta_k\) parret med sig selv, dvs. integralet

\[\begin{equation*} >K^{(k)}_{11} = \int_{x_k}^{x_{k+1}} EI(x)\,\zeta_k''(x)^2\,dx. >\end{equation*}\]

Beregn dette integral for konstant \(EI = 1\) og eftervis, at

\[\begin{equation*} K^{(k)}_{11} = \frac{12\,EI}{h^3}. \end{equation*}\]

Opgave D.5#

Implementer samlingen af den globale \(2N \times 2N\) stivhedsmatrix for bjælkeligningen og løs systemet for \(N = 6\) med \(EI = 1\) og \(P = 1\). Plot FEM-løsningen sammen med den eksakte løsning fra Opgave D.1.

Opgave D.6#

Løs bjælkeligningen for \(N = 2, 4, 8, 16\) og vis i fire delplot FEM-løsningen sammen med den eksakte løsning. Vurder på den baggrund, om løsningen synes at konvergere hurtigere eller langsommere end for andenordens-problemet.

Udvidelse E - Ikke-homogene randbetingelser#

Randbetingelserne \(u(0) = u(1) = 0\), som vi har brugt gennem hele projektet, er et specielt tilfælde. I mange fysiske situationer er randværdierne ikke nul, for eksempel en stang hvis endepunkter holdes på forskellige temperaturer, eller en snor der er fastgjort i forskellige højder. Matematisk skrives dette som \(u(0) = \alpha\) og \(u(1) = \beta\) med \(\alpha, \beta\) ikke nødvendigvis lig nul. Spørgsmålet er, hvordan FEM-metoden tilpasses.

Standardtricket er en løftningsfunktion: Skriv

\[ u(x)=\tilde{u}(x)+L(x), \]

hvor \(L(x)=\alpha+(\beta-\alpha)x\) er den lineære funktion, der opfylder randbetingelserne, og hvor \(\tilde{u}(x)\) er den ubekendte korrektion med

\[ \tilde{u}(0)=\tilde{u}(1)=0. \]

Opgave E.1#

(a) Udled det modificerede lineære system for \(\tilde{u}\) og identificer den ændrede lastvektor.

(b) Løs numerisk for \(a(x) = 1\), \(p(x) = 1\), \(\alpha = 1\), \(\beta = 0\) med \(N = 20\). Plot løsningen og sammenlign med den eksakte løsning, som kan findes ved direkte integration.

Contents