Elementtimenetelmän kantafunktioiden laskeminen

Esitellään toimiva tapa rakentaa elementin kantafunktiot.

Aiemmin ei juurikaan menty siihen, kuinka päädyttiin esitettyihin kantafunktioihin \[ \begin{alignedat}{1} \phi_{2}^{\left(1\right)} & =\frac{1}{2}x,\\ \phi_{1}^{\left(2\right)} & =3-x. \end{alignedat} \]

Tässä tapauksessa tehtävä on vielä niin helppo, että sen voi jopa päätellä. Katsotaan nyt kuitenkin hieman systemaattisempi tapa hoitaa homma. Lähtökohtana on muodostaa elementin alueelle interpolaatiopolynomit \(\phi_i(x)\) sopivien ehtojen mukaisesti. Klassisen interpolaatiopolynomien perheen muodostaa Lagrangen elementit, jotka saavat arvon 1 kyseissä vapausasteessa \(x_i\) ja arvon 0 kaikissa muissa vapausasteissa.

Ensimmäisessä lineaarisessa \(\mathcal{C}^0\) elementissä laskenta-alue olisi välillä \(\Omega^{(1)} = [0, 2]\). Polynomit ovat joukossa \({\{1, x\}}\). Vapausaste \(x_1 = 0\) ja vapausaste \(x_2\) = 2. Tästä voidaan rakentaa Vandermonden matriisi \[ \boldsymbol{V}_{1}=\left[\begin{array}{cc} 1 & 0\\ 1 & 2 \end{array}\right], \]

Jonka käänteismatriisista saadaan kertoimet \[ \begin{alignedat}{1} \boldsymbol{c}_{1}^{\left(1\right)} & =\boldsymbol{V}_{1}^{-1}\left[\begin{array}{c} 1\\ 0 \end{array}\right]=\left[\begin{array}{r} 1\\ -\frac{1}{2} \end{array}\right]\\ \boldsymbol{c}_{2}^{\left(2\right)} & =\boldsymbol{V}_{1}^{-1}\left[\begin{array}{c} 0\\ 1 \end{array}\right]=\left[\begin{array}{c} 0\\ \frac{1}{2} \end{array}\right] \end{alignedat} \]

Elementin \(\Omega^{(1)}\) muotofunktiot siis ovat \[ \begin{alignedat}{1} \phi_{1}^{\left(1\right)}\left(x\right) & =1-\frac{1}{2} x\\ \phi_{2}^{\left(1\right)}\left(x\right) & =\frac{1}{2}x \end{alignedat} \]

Vastaavasti toiselle elementille saadaan Vandermonden matriisi \[ \begin{alignedat}{1} \boldsymbol{V}_{2}=\left[\begin{array}{cc} 1 & 2\\ 1 & 3 \end{array}\right] \end{alignedat} \]

Elementin \(\Omega^{(2)}\) muotofunktiot ovat siten \[ \begin{alignedat}{1} \phi_{1}^{\left(2\right)}\left(x\right) & =3-x\\ \phi_{2}^{\left(2\right)}\left(x\right) & =x-2 \end{alignedat} \]

Muotofunktiot voidaan laskea symbolisesti esimerkiksi SymPy:llä:

from sympy import *
var('x')

P1 = Matrix([1, x])
X1 = Matrix([0, 2])

V1 = zeros(2, 2)
for i in range(2):
    for j in range(2):
        V1[i,j] = P[j].subs({x: X1[i]})

c1 = V1.inv() * Matrix([1, 0])
c2 = V1.inv() * Matrix([0, 1])
phi1_1 = (c1.T * P1)[0]
phi1_2 = (c2.T * P2)[0]

print(phi1_1)
print(phi1_2)
1 - x/2
x/2

Saimme tässä prosessissa oikeastaan neljä muotofunktiota, jotka liittyvät kahteen elementtiin. Kuitenkin kiinnostuksen kohteena oli vain keskimmäinen solmupiste, sillä tehtävässä oli homogeeniset reunaehdot ja siten ensimmäiseen ja viimeiseen solmupisteeseen liittyvät muotofunktiot voitiin poistaa ongelmasta.

Samanlaisella tekniikalla saamme rakennettua esimerkiksi kuubiset \(\mathcal{C}^1\) jatkuvat Hermiten interpolaatiopolynomit. Nyt polynomien joukko on \[ P=\left\{ 1,x,x^{2},x^{3}\right\} \] Ja vapausasteita tarvitsemme neljä: funktion ja sen derivaatan jatkuvuudet elementin molemmissa päissä. Vandermonden matriisin rakentaminen ja interpolaatiopolynomit ovat

P1 = Matrix([1, x, x**2, x**3])
dP1 = P1.diff(x)
X1 = Matrix([0, 0, 2, 2])

V1 = Matrix([
    P1.subs({x: X1[0]}).T,
    dP1.subs({x: X1[1]}).T,
    P1.subs({x: X1[2]}).T,
    dP1.subs({x: X1[3]}).T
])

phi1_1 = ((V1.inv() * Matrix([1, 0, 0, 0])).T * P1)[0]
phi2_1 = ((V1.inv() * Matrix([0, 1, 0, 0])).T * P1)[0]
phi3_1 = ((V1.inv() * Matrix([0, 0, 1, 0])).T * P1)[0]
phi4_1 = ((V1.inv() * Matrix([0, 0, 0, 1])).T * P1)[0]

\[ \begin{alignedat}{1} \phi_{1}^{\left(1\right)}\left(x\right) & =\frac{1}{4}x^{3}-\frac{3}{4}x^{2}+1\\ \phi_{2}^{\left(1\right)}\left(x\right) & =\frac{1}{4}x^{3}-x^{2}+x\\ \phi_{3}^{\left(1\right)}\left(x\right) & =-\frac{1}{4}x^{3}+\frac{3}{4}x^{2}\\ \phi_{4}^{\left(1\right)}\left(x\right) & =\frac{1}{4}x^{3}-\frac{1}{2}x^{2} \end{alignedat} \]

Vastaavasti alueelle \(\Omega^{(2)} = [2, 3]\) saadaan

\[ \begin{alignedat}{1} \phi_{1}^{\left(2\right)}\left(x\right) & =2x^{3}-15x^{2}+36x-27\\ \phi_{2}^{\left(2\right)}\left(x\right) & =x^{3}-8x^{2}+21x-18\\ \phi_{3}^{\left(2\right)}\left(x\right) & =-2x^{3}+15x^{2}-36x+28\\ \phi_{4}^{\left(2\right)}\left(x\right) & =x^{3}-7x^{2}+16x-12 \end{alignedat} \]

Kaikki muotofunktiot eivät ole kinemaattisesti käypiä. Kinemaattisesti käyvät on esitetty kuvassa 1.

Figure 1: Differentiaaliyhtälön tarkka ratkaisu.

Sijoittamalla heikkoon muotoon saadaan

\[ \boldsymbol{K}=\left[\begin{array}{cc} \frac{27}{2} & \frac{9}{2}\\ \frac{9}{2} & 6 \end{array}\right] \]

\[ \boldsymbol{f}=\left[\begin{array}{r} 36\\ -6 \end{array}\right] \]

\[ \boldsymbol{u}=\left[\begin{array}{c} u_{2}\\ v_{2} \end{array}\right]=\boldsymbol{K}^{-1}\boldsymbol{f}=\left[\begin{array}{r} 4\\ -4 \end{array}\right] \]

Ratkaisu on esitetty kuvassa 2.

Figure 2: Differentiaaliyhtälön ratkaisu elementtimenetelmällä ja C1-jatkuvilla Hermiten kuubisilla muotofunktioilla.