MES równania różniczkowe

Przykład

Wyznaczyć rozwiązanie MES dyskretyzując dziedzinę jednym elementem skończonym z kwadratowymi hierarchicznymi funkcjami kształtu dla problemu brzegowego:

\begin{aligned} &y''(x) - 2y + 2x(x - 3) = 0, \quad x \in [-1, 2] \\ &y(-1) = 5 \\ &y(2) = -1 \end{aligned}

1. Obliczenie resisduum

Dla podanej postaci jest to po prostu:

\begin{aligned} &R(x)=y''(x) - 2y + 2x(x - 3) \\ \end{aligned}

2. Przemnożenie residuum przez funkcję wagową

\begin{aligned} &R(x)\cdot \omega(x)=y''(x)\cdot \omega(x) - 2y\cdot \omega(x) + 2x(x - 3)\cdot \omega(x) \\ \end{aligned}

3. Minimalizacja residuum

Tutaj należy się na chwilę zatrzymać, we wstępie granice całkowania były podane jako \(x_a, x_b\), natomiast my będziemy każdy element zaczynali w \(x_a=0\) korzystając z faktu że możemy układać zaczepić lokalnie w punkcie początkowym danego elementu skończonego.

\begin{aligned} \int_{0}^{l} \omega \cdot y''(x) \, dx - \int_{0}^{l} \omega \cdot 2y \, dx + \int_{0}^{l} \omega \cdot 2x \cdot (x - 3) \, dx = 0 \end{aligned}

Musimy przekształcić pierwszy element wyrażenia, skorzystamy z całkowania przez części:

\begin{aligned} &\int u \, dv = u \cdot v - \int v \, du \\ &\int_{0}^{l} \omega \cdot y'' \, dx = \omega \cdot y' \Big{|}_{0}^{l} - \int_{0}^{l} \omega' \cdot y' \, dx \end{aligned}

Wracamy do równania:

\begin{aligned} \omega \cdot y' \Big{|}_{0}^{l} - \int_{0}^{l} \omega' \cdot y' \, dx - \int_{0}^{l} \omega \cdot 2y \, dx + \int_{0}^{l} \omega \cdot 2x \cdot (x - 3) \, dx = 0 \end{aligned}

4. Aproksymacja

Korzystając z poprzednio podanych zależności przekształcamy równanie:

\begin{aligned} y_e &= N_e \cdot d \quad &y'_e &= N'_e \cdot d \\ \omega_e &= \beta_e^{T} \cdot N_e^{T} \quad &\omega'_e &= \beta_e^{T} \cdot N{'}_e^{T} \end{aligned} \begin{aligned} &\int_{0}^{l} \omega' y' \, dx + \int_{0}^{l} \omega \cdot 2y \, dx = \omega \cdot y' \Big|_{0}^{l} + \int_{0}^{l} \omega \cdot 2x \cdot (x - 3) \, dx \\ &\int_{0}^{l} \beta_e^T N_e^T N'_e \cdot d \, dx + \int_{0}^{l} \beta_e^T N_e^T \cdot 2 N_e \cdot d \, dx = \beta_e^T N_e^T y' \Big|_{0}^{l} + \int_{0}^{l} \beta_e^T N_e^T \cdot 2x \cdot (x - 3) \, dx \\ &\int_{0}^{l} N{'}_e^T N'_e \cdot d \, dx + \int_{0}^{l} N_e^T \cdot 2 N_e \cdot d \, dx = N_e^T y' \Big|_{0}^{l} + \int_{0}^{l} N_e^T \cdot 2x \cdot (x - 3) \, dx \end{aligned}

Przy okazji podzieliliśmy równanie na lewą i prawą stronę i uzyskaliśmy postać którą można zapisać jako:

\begin{aligned} K\cdot d=f_b+f \end{aligned}

Mamy jeden element skończony, zapisujemy dla niego:

początek elementu w układzie globalnym \(\quad a = -1\)

długość elementu \(\quad l = 3\)

\begin{aligned} x_e &= x - a \quad \Rightarrow \quad x_e = x + 1 \\ x &= x_e + a \quad \Rightarrow \quad x = x_e - 1 \end{aligned}

Obliczamy poszczególne składniki równania. Przyjęta aproksymacja daje nam:

\begin{aligned} N_e(x_e) &= \begin{bmatrix} 1 - \frac{x_e}{l} & \frac{x_e}{l} \end{bmatrix} x_e \cdot (x_e - l) \quad \Rightarrow \quad \begin{bmatrix} -\frac{x_e}{3} + 1 & \frac{x_e}{3} x_e \cdot (x_e - 3) \end{bmatrix} \\ N'_e(x_e) &= \frac{d}{dx_e} N_e(x_e) \quad \Rightarrow \quad \begin{bmatrix} -\frac{1}{3} & \frac{1}{3} \cdot 2 \cdot x_e - 3 \end{bmatrix} \end{aligned}

Więc lewa strona:

\begin{aligned} \mathbf{K} &= \int_{0}^{l} N'_e(x_e) ^T \cdot N'_e(x_e) \, dx_e + \int_{0}^{l} N_e(x_e) ^T \cdot 2 \, N_e(x_e) \, dx_e \\ \end{aligned} \begin{aligned} &K=\int_0^l\left[\begin{array}{c} -\frac{1}{3} \\ \frac{1}{3} \\ 2 \cdot x_e-3 \end{array}\right] \cdot\left[-\frac{1}{3} \frac{1}{3} 2 \cdot x_e-3\right] \mathrm{d} x_e+\int_0^l\left[\begin{array}{c} -\frac{x_e}{3}+1 \\ \frac{x_e}{3} \\ x_e \cdot\left(x_e-3\right) \end{array}\right] \cdot 2\left[-\frac{x_e}{3}+1 \frac{x_e}{3} x_e \cdot\left(x_e-3\right)\right] \mathrm{d} x_e \rightarrow\left[\begin{array}{ccc} \frac{7}{3} & \frac{2}{3} & -\frac{9}{2} \\ \frac{2}{3} & \frac{7}{3} & -\frac{9}{2} \\ -\frac{9}{2} & -\frac{9}{2} & \frac{126}{5} \end{array}\right] \end{aligned}

Oraz d zgodnie z warunkami zadania:

\begin{aligned} d &= \begin{bmatrix} d_1 \\ d_2 \\ \alpha \end{bmatrix} \quad \Rightarrow \quad \begin{bmatrix} 5 \\ -1 \\ \alpha \end{bmatrix} \end{aligned}

Prawa strona, uwaga - pamiętamy o zamianie \(x\) globalnego na \(x_e\) lokalne:

\begin{aligned} p &= \int_{0}^{l} \left( N_e(x_e) \right)^T \cdot 2x \cdot (x - 3) \, dx_e \\ p &= \int_{0}^{l} \begin{bmatrix} \frac{x_e}{3} + 1 \\ \frac{x_e}{3} \\ x_e \cdot (x_e-3) \\ \end{bmatrix} \cdot 2 \cdot (x_e + a) \cdot (x_e + a - 3) \, dx_e \quad \Rightarrow \quad \begin{bmatrix} \frac{3}{2} \\ \frac{9}{2} \\ \frac{36}{5} \end{bmatrix} \end{aligned}

Drugi element prawej strony jest trywialny:

\begin{aligned} p_b = \begin{bmatrix} -y'(0) \\ y'(l) \\ 0 \end{bmatrix} \end{aligned}

Ostatecznie:

\begin{aligned} \mathbf{K} \cdot d &= p_b + p \\ \begin{bmatrix} \frac{7}{3} & \frac{2}{3} & \frac{9}{2} \\ \frac{2}{3} & \frac{7}{3} & \frac{9}{2} \\ \frac{9}{2} & \frac{9}{2} & \frac{126}{5} \end{bmatrix} \cdot \begin{bmatrix} 5 \\ -1 \\ \alpha \end{bmatrix} &= \begin{bmatrix} -y'(0) \\ y'(l) \\ 0 \end{bmatrix} + \begin{bmatrix} \frac{3}{2} \\ \frac{9}{2} \\ \frac{36}{5} \end{bmatrix} \end{aligned}

Rozwiązanie układu daje nam:

\begin{aligned} y'(0) &= -5 \\ y'(l) &= 1 \\ \alpha &= 1 \end{aligned}

Stąd d:

\begin{aligned} d &= \begin{bmatrix} d_1 \\ d_2 \\ \alpha \end{bmatrix} \quad \Rightarrow \quad \begin{bmatrix} 5 \\ -1 \\ 1 \end{bmatrix} \end{aligned}

Więc rozwiązanie przybliżone w układzie lokalnym i globalnym:

\begin{aligned} y_p &= N_e(x_e) \cdot d \\ y_p &= \begin{bmatrix} 1 - \frac{x_e}{l} & \frac{x_e}{l} \end{bmatrix} x_e \cdot (x_e - l) \cdot \begin{bmatrix} 5 \\ -1 \\ 1 \end{bmatrix} = x_e^2 - 5 \cdot x_e + 5 \\ x_e &\equiv x - a \\ y_p &= x_e^2 - 5 \cdot x_e + 5 = x^2 - 3 \cdot x + 1 \end{aligned}