Странный результат при решении МКЭ

[Geckelberryfinn]

Стоит задача решить уравнение (остро,блин стоит-- диплом таки)
\mathcal{L} u=\chi(x,t),
В симметричной области, ограниченной двумя окружностями (т.е. в кольце).
Оператор уравнения и ГУ выглядят так:
\mathcal{L}=\left(x^2\frac{\partial^2}{\partial x^2}+y^2\frac{\partial^2}{\partial y^2}+xy\frac{\partial}{\partial x}\frac{\partial}{\partial y}+x\frac{\partial}{\partial x} + y\frac{\partial}{\partial y} -1\right)}\frac{1}{(x^2+y^2)}
\left((\kappa+1)\left[x\frac{\partialu}{\partial x}+y\frac{\partial u}{\partial y}\right]+(\kappa-1)u\right)\frac{1}{\sqrt{x^2+y^2}}=q(x,y,t), если x \in S_1 (внешняя граница)
u=0, если x \in S_2 (внутри)
Решаю на треугольной сетке согласно Зинкевичу, минимизируя невязку по области и на границе
\int\limits_{\Omega} N_i R_\Omega d\Omega+\int\limits_{\Gamma} N_i R_\Gamma d\Gamma=0
где N_i - базисные, они же пробные функции, а R_\Omega, R_\Gamma - невязки по области и по границе соотвественно.
Тогда элементы матрицы левой части СЛАУ и ветор правой части, я вычисляю следующим образом
K_{i,j}=\int\limits_{\Omega} N_i\mathcal{L} N_j d\Omega+\int\limits_{\Gamma} N_i\mathcal{M}N_j d\Gamma, f_i=-\int\limits_{\Omega} N_i \chi d\Omega-\int\limits_{\Gamma} rd\Gamma
где M и r таковы, что на границе
\mathcal{M}u=r
Перед тем, как интегрировать, интеграл по области, записанный в выражении для элементов матрицы левой части, взял по частям, для того, чтобы избавиться от вторых производных (посколько базисные функции должны в таком случае принадлежать не C^1, а C^0.
Интегрирование по треугольникам производится пока методом Гаусса по семи точкам. По поводу этого алгоритма могу сказать, что, скорее всего, он верен. Советовали аналитически интеграл вычислять, но пока боюсь. Вроде бы хоть и маленькая, но уверенность в правильной работе алгоритма интегрирования есть.
Далее, базисные функции - квадратичные. Внешний вид прилагаю. Вроде бы тоже верны.
Сетка тоже вроде хорошая. При ее сгущении качественно ничего не меняется. СЛАУ решается LU-алгоритмом. (также еще и QR пробовалось. Вроде тоже верно, с MathCad сходится)
Много вопросов возникает при виде получающегося решения.
Может быть кто-нибудь сталкивался с таким поведением. Из-за чего это может происходить?
Буду очень признателен за любую помощь.
Базисные функции


Сетка

Решение

[Zai]

Я не помню точно требования к дипломным работам, но первое что нужно сделать это подготовить оператор к требованиям МКЭ

Нужно домножить на (x^2+y^2) . Уравнение у Вас что-то вроде Лапласа с внутренним источником тепла и переносом. Не нужно перемалывать функциями формы вектор внешних сил - вдруг будет ошибка. Помножте среднее значение на площадь треугольника и разделите на 3 - вот вам и вектор правых частей для тругольника.


Граничное условие у Вас содержит градиент функции по нормали к границе и тоже достаточно просто(если Вы перейдете в локальную систему координат то X и Y исчезнут).

В МКЭ методах основная гипотеза - постоянство коэффициентов оператора внутри треугольника. Если размеры элементов достаточно малы то переменность коэффициентов также аппроксимируется. Как Вы это используете?

Есть также правило. Если у Вас новый оператор, то перед тем как его решать, решите сначала с простым оператором Лапласа с простыми граничными условиями и сравните с аналитическим решением,

То же самое с граничными условиями. Задайте на внешней границе 1 и убедитесь что решение физично.

[Zai]

На практике если Вам необходимо получть решение за короткий срок, лучше воспользоваться линейными функциями формы.
Минимизируемый по пространству функционал будет иметь вид:
$2W=X^2(\frac {\partial u} {\partial x})^2+Y^2(\frac {\partial u} {\partial y})^2+XY\frac {\partial u} {\partial x}\frac {\partial u} {\partial y}+Xu\frac {\partial u} {\partial x}+Yu\frac {\partial u} {\partial y}- {u^2}
Отмечаю еще раз что коэффициенты постоянны внутри каждого треугольника.
При линейных функциях формы сразу получаем, что

$ u=\frac {u_1+u_2+u_3} 3
$\frac {\partial u} {\partial x}=\alpha_1u_1+\alpha_2u_3+\alpha_3u_3
$\frac {\partial u} {\partial y}=\beta_1u_1+\beta_2u_3+\beta_3u_3

Значения производных получены из первых коэффициентов разложения в ряд Тейлора:
$ u_2=u_1+\frac {\partial u} {\partial x}*(x_2-x_1)+\frac {\partial u} {\partial y}*(y_2-y_1)
$ u_3=u_1+\frac {\partial u} {\partial x}*(x_3-x_1)+\frac {\partial u} {\partial y}*(y_3-y_1)
Перемножая W*S и дифференцииируя выражение по $u_1 ,u_2 ,u_3 Вы получите матрицу жесткости.
Изложенная выше процедура может быть проведена вручную и Вы можите проверить матрицу жесткости одного из Ваших элементов. Они не должны значительно отличаться для достаточно мелких элементов.
С граничными условиями на второй границе сложнее. Так как когда решается уравнение математичекой физики мы можем однозначно связать энерию тела и работу внешних сил.
У Вас на ребре
$u=\frac {u_2+u_3} 2 если ребро 2-3 граничное.
Приставим к ребру квадрат со сторонами ребра 2-3
и запишем минимизируемое выражение
$W_1l^2= \left(\left((\kappa+1)\left[x\frac{\partial u}{\partial x}+y\frac{\partial u}{\partial y}\right]+(\kappa-1)u\right)\frac{1}{\sqrt{x^2+y^2}}-q(x,y,t)\right)^2
Дифференциируя по $u_1 ,u_2 ,u_3 получим добавочные коэффициенты как для матрицы жесткости так и добавочные значения правых частей.