Простейший пример МКЭ (Галеркин)

[_genius_]

https://homepages.cae.wisc.edu/~suresh/ ... trofem.pdf учебное пособие, самое начало 1.3.1
https://www.twirpx.com/file/627919/ книга от того же автора, эта глава в ней один к одному как в пособии только под номером 1.2.1

Есть одномерное ДУ 2го порядка:

$a \frac{d^2 u}{dx^2} + b = 0, \qquad 0 \leq x \leq 2L$

с граничными условиями:

$\begin{cases} u \vert_{x=0} = 0 \\ a\frac{du}{dx} \vert_{x=2L} = R \end{cases}$

Его очень легко решить аналитически, что я собственно и сделал для того чтобы вручную проверить каждый шаг:

$u = \frac{-b}{2a}x^2 + (\frac{2bL}{a} + \frac{R}{a})x$

$\frac{du}{dx} = \frac{-b}{a}x + (\frac{2bL}{a} + \frac{R}{a})$

Область определения делится на два равных участка (от 0 до L и от L до 2L) и на каждом вводятся по две линейные базисные функции ("shape function") через которые с помощью узловых значений приближается значение искомой функции u на данном участке.

Для первого участка:

$u_{12} = N_1 u_1 + N_2 u_2$

$N_1 = \frac{x_2-x}{x_2-x_1}; \qquad N_2 = \frac{x-x_1}{x_2-x_1}$

Затем эта аппроксимация подставляется в оригинальное уравнение, оно домножается на вектор с базисными функциями и интегрируется (невязка, условие ортогональности и т.д., это я все понимаю):

$\int_{x_1}^{x_2} \mathbf{N}^T a \frac{d^2}{dx^2} \mathbf{N} \mathbf{u} dx + \int_{x_1}^{x_2} \mathbf{N}^T b dx = 0$

Векторы $\mathbf{u} = \begin{bmatrix} u_1 \\ u_2 \end{bmatrix}$ и $\mathbf{N} = \begin{bmatrix} N_1 & N_2 \end{bmatrix}$
А вот тут становится интересно, автор с помощью интегрирования по частям приходит к такому результату:

$\int_{x_1}^{x_2} \frac{d}{dx} \mathbf{N}^T a \frac{d}{dx} \mathbf{N} dx \mathbf{u} - \int_{x_1}^{x_2} \mathbf{N}^T b dx - \begin{bmatrix} 0 \\ 1 \end{bmatrix} a \frac{du}{dx}\vert_{x=x_2} + \begin{bmatrix} 1 \\ 0 \end{bmatrix} a \frac{du}{dx}\vert_{x=x_1} = 0 $

Первый член обозначается как $k u$, остальные как $f$ и переносятся вправо, затем подставляются пределы и для первого участка система выглядит так:

$\frac{a}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \frac{bL}{2} \begin{bmatrix} 1 \\ 1 \end{bmatrix}$

А для второго участка так:

$\frac{a}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} = \frac{bL}{2} \begin{bmatrix} 1 \\ 1 \end{bmatrix} + \begin{bmatrix} 0 \\ R \end{bmatrix}$

Что тут странного: при вычислении k вектор с узловыми точками был фактически вынесен за скобку, что в общем-то правильно, т.к. от x он не зависит, однако при вычислении составной части f приближенная функция $\mathbf{N} \mathbf{u}$ под дифференциалом была заменена обратно на оригинальную функцию $u$. Этот трюк позволил в последнем случае получить вектор, а не матрицу. И все бы ничего, но значение производной известно только в третьей узловой точке, но не в первых двух. Судя по итоговому значению f на первом участке автор просто взял и приравнял их к нулю. Одного взгляда на аналитическое решение достаточно чтобы понять что это неправильно.

Далее процедура повторяется для второго участка, k и f накладываются друг на друга во втором узле, получается СЛАУ из 3х уравнений, затем применяется первое граничное условие чтобы обнулить коэффициенты в первой строчке:

$\frac{a}{L} \begin{bmatrix} 1 & 0 & 0\\ 0 & 2 & -1 \\ 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \frac{bL}{2} \begin{bmatrix} 0 \\ 2 \\ 1 \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ R \end{bmatrix}$

Эта система абсолютно верная и приводит к точным значениям u в трех узлах:

$u_2 = \frac{1.5b L^2}{a} + \frac{r L}{a}$

$u_3 = \frac{2b L^2}{a} + \frac{2r L}{a}$

Разумеется мне захотелось узнать как это могло получиться и поэтому я перемножил k и $\mathbf{u}$ вручную на обоих участках чтобы получить соответствующие f, в результате:

$\frac{L}{a} f_1 = \begin{bmatrix} -\frac{1.5b L^2}{a} - \frac{r L}{a} \\ \frac{1.5b L^2}{a} + \frac{r L}{a} \end{bmatrix}$

$\frac{L}{a} f_2 = \begin{bmatrix} -\frac{0.5b L^2}{a} - \frac{r L}{a} \\ \frac{0.5b L^2}{a} + \frac{r L}{a} \end{bmatrix}$

Очевидно что при применении первого граничного условия первая строчка $f_1$ исчезнет, так что она не имеет значения, к последней строчке вопросов нет, она совпадает с авторской.
А вот со второй строчкой $f_1$ и первой строчкой $f_2$ все куда интереснее - слагаемое содержащее R исчезает только после их сложения.

И в самом деле если подставить значения производной в первых двух узлах $u'(x_1) = \frac{2bL}{a} + \frac{R}{a}$ и $u'(x_2) = \frac{bL}{a} + \frac{R}{a}$, которые автор беспардонно обнулил в формулу для f:

$\int_{x_1}^{x_2} \mathbf{N}^T b dx + \begin{bmatrix} 0 \\ 1 \end{bmatrix} a \frac{du}{dx}\vert_{x=x_2} - \begin{bmatrix} 1 \\ 0 \end{bmatrix} a \frac{du}{dx}\vert_{x=x_1} = \frac{bL}{2} \begin{bmatrix} 1 \\ 1 \end{bmatrix} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} a \frac{du}{dx}\vert_{x=x_2} - \begin{bmatrix} 1 \\ 0 \end{bmatrix} a \frac{du}{dx}\vert_{x=x_1}$

получаются как раз вышенаписанные $f_1$ и $f_2$.

Это что же получается, откровенный мухлеж от профессора с кучей цитирований https://scholar.google.com/citations?user=SxTP4dMAAAAJ который не заметили ни студенты на лекциях в Японии, ни редакторы книги?

У меня уже 5ый день крыша едет от всего этого, уже десятки раз все перепроверил. Можно ли как-то подправить этот пример чтобы честно прийти к той же системе или это гиблое дело? Почему во всем интернете нет ни одного легкого примера по МКЭ который можно было бы объяснить любому идиоту на пальцах? Либо как здесь с ошибками и мухлежом, либо такая жесть, которую поймут только доктора физ.мат. наук, да и то не все. Главное работать с пакетами очень просто, но вот объяснить что же там на самом деле происходит и почему это так практически невозможно.

[_genius_]

В общем помог себе сам: с такими базисными функциями решить задачу без мухляжа невозможно. Их необходимо заменить на те, которые удовлетворяют однородному граничному условию, т.е. $N_1(0) = N_2(0) = 0$, к примеру $N_1 = x, N_2 = x^2$. Более того, двух базисных функций достаточно, придумывать по две для каждого элемента не нужно.

[_genius_]

Но правда при этом получается что элементов вроде бы и нет. Чтобы они были можно задать кусочные базисные функции, но опять-таки соблюдая граничное условие, тогда тоже все получается.
Вот $N_1$:
$\begin{cases} \frac{x}{L}, 0 \le x \le L \\ \frac{2L-x}{L}, L \le x \le 2L\end{cases}$
Вот $N_2$:
$\begin{cases} 0, 0 \le x \le L \\ \frac{x-L}{L}, L \le x \le 2L\end{cases}$