J. Пример применения функционала

Рассмотрим, как функционал метода Бубнова–Галёркина из приложения I применяется на практике: разберём одномерный пример от начала до конца — от пробных функций до системы линейных уравнений в матричной форме.

Решаем уравнение теплопроводности L[u]L[u]​=f{} = f, где параболический оператор LL​=∂u∂t{} = \frac{\partial u}{\partial t}​−a2Δu{} - a^2 \Delta u, на отрезке [x0,x4][x_0, x_4]. На концах зададим граничные условия Неймана (1.14) — заданы потоки ∂u∂n∣x=x0\frac{\partial u}{\partial n}\bigg|_{x=x_0}​=g0{} = g_0 и ∂u∂n∣x=x4\frac{\partial u}{\partial n}\bigg|_{x=x_4}​=g4{} = g_4. Доведём функционал Бубнова–Галёркина до матричной формы и подставим эти условия.

Одномерная сетка из четырёх симплексов-отрезков s₁…s₄ и пяти узлов
Рис. J.1. Одномерная сетка примера: четыре симплекса-отрезка s1,…,s4s_1, \ldots, s_4 и пять узлов x0,…,x4x_0, \ldots, x_4; длины отрезков обозначаются l(i)(i+1)l_{(i)(i+1)}.

Возьмём четыре симплекса-отрезка s1,s2,s3,s4s_1, s_2, s_3, s_4, три внутренние точки x1,x2,x3x_1, x_2, x_3 и две граничные точки x0,x4x_0, x_4 (рисунок J.1). Пробные функции для симплексов-отрезков примут вид

υ(0)(1)(x)\displaystyle \upsilon_{(0)(1)}(x) =q0⋅ϕ0(x)\displaystyle {} = q_0 \cdot \phi_0(x)​+q1⋅ϕ1(x)\displaystyle {} + q_1 \cdot \phi_1(x)υ(1)(2)(x)\displaystyle \upsilon_{(1)(2)}(x) =q1⋅ϕ1(x)\displaystyle {} = q_1 \cdot \phi_1(x)​+q2⋅ϕ2(x)\displaystyle {} + q_2 \cdot \phi_2(x)υ(2)(3)(x)\displaystyle \upsilon_{(2)(3)}(x) =q2⋅ϕ2(x)\displaystyle {} = q_2 \cdot \phi_2(x)​+q3⋅ϕ3(x)\displaystyle {} + q_3 \cdot \phi_3(x)υ(3)(4)(x)\displaystyle \upsilon_{(3)(4)}(x) =q3⋅ϕ3(x)\displaystyle {} = q_3 \cdot \phi_3(x)​+q4⋅ϕ4(x)\displaystyle {} + q_4 \cdot \phi_4(x)

Функции «крышки» и их производные для симплексов-отрезков имеют вид

ϕ0(x)\displaystyle \phi_0(x) =a0\displaystyle {} = a_0​+b0⋅xϕ0′(x)\displaystyle {} + b_0 \cdot x \quad \phi_0'(x)​=b0\displaystyle {} = b_0ϕ1(x)\displaystyle \phi_1(x) =a1\displaystyle {} = a_1​+b1⋅xϕ1′(x)\displaystyle {} + b_1 \cdot x \quad \phi_1'(x)​=b1\displaystyle {} = b_1ϕ2(x)\displaystyle \phi_2(x) =a2\displaystyle {} = a_2​+b2⋅xϕ2′(x)\displaystyle {} + b_2 \cdot x \quad \phi_2'(x)​=b2\displaystyle {} = b_2ϕ3(x)\displaystyle \phi_3(x) =a3\displaystyle {} = a_3​+b3⋅xϕ3′(x)\displaystyle {} + b_3 \cdot x \quad \phi_3'(x)​=b3\displaystyle {} = b_3ϕ4(x)\displaystyle \phi_4(x) =a4\displaystyle {} = a_4​+b4⋅xϕ4′(x)\displaystyle {} + b_4 \cdot x \quad \phi_4'(x)​=b4\displaystyle {} = b_4

Коэффициенты принимают следующие значения. Введём обозначение l(i)(i+1)l_{(i)(i+1)}​=xi+1{} = x_{i+1}​−xi{} - x_i для длины симплекса между точками xix_i и xi+1x_{i+1}.

Коэффициенты «крышек» для произвольного симплекса между точками xix_i и xi+1x_{i+1} выводятся так же, как в разделе «Пробные функции», и равны

ai\displaystyle a_i =xi+1l(i)(i+1)bi\displaystyle {} = \frac{\displaystyle x_{i+1}}{\displaystyle l_{(i)(i+1)}} \quad b_i​=−1l(i)(i+1)\displaystyle {} = \frac{\displaystyle -1}{\displaystyle l_{(i)(i+1)}}ai+1\displaystyle a_{i+1} =−xil(i)(i+1)bi+1\displaystyle {} = \frac{\displaystyle -x_i}{\displaystyle l_{(i)(i+1)}} \quad b_{i+1}​=1l(i)(i+1)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(i)(i+1)}}

Для краевых «крышек» ϕ0\phi_0 и ϕ4\phi_4 берётся соответствующая половина. Подставив эти коэффициенты, сразу запишем пробные функции на каждом симплексе.

Подставим коэффициенты в пробные функции и выразим их относительно xx:

Для симплекса s1s_1:

υ(0)(1)(x)\displaystyle \upsilon_{(0)(1)}(x) =q0⋅ϕ0(x)\displaystyle {} = q_0 \cdot \phi_0(x)​+q1⋅ϕ1(x)\displaystyle {} + q_1 \cdot \phi_1(x)​=q0\displaystyle {} = q_0​⋅(x1l(0)(1)−xl(0)(1))\displaystyle {} \cdot \left(\frac{\displaystyle x_1}{\displaystyle l_{(0)(1)}} - \frac{\displaystyle x}{\displaystyle l_{(0)(1)}}\right)​+q1\displaystyle {} + q_1​⋅(−x0l(0)(1)+xl(0)(1))\displaystyle {} \cdot \left(\frac{\displaystyle -x_0}{\displaystyle l_{(0)(1)}} + \frac{\displaystyle x}{\displaystyle l_{(0)(1)}}\right)​=1l(0)(1)[q0⋅(x1−x)+q1⋅(x−x0)]\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} \left[q_0 \cdot (x_1 - x) + q_1 \cdot (x - x_0)\right]

Для симплекса s2s_2:

υ(1)(2)(x)\displaystyle \upsilon_{(1)(2)}(x) =q1⋅ϕ1(x)\displaystyle {} = q_1 \cdot \phi_1(x)​+q2⋅ϕ2(x)\displaystyle {} + q_2 \cdot \phi_2(x)​=q1\displaystyle {} = q_1​⋅(x2l(1)(2)−xl(1)(2))\displaystyle {} \cdot \left(\frac{\displaystyle x_2}{\displaystyle l_{(1)(2)}} - \frac{\displaystyle x}{\displaystyle l_{(1)(2)}}\right)​+q2\displaystyle {} + q_2​⋅(−x1l(1)(2)+xl(1)(2))\displaystyle {} \cdot \left(\frac{\displaystyle -x_1}{\displaystyle l_{(1)(2)}} + \frac{\displaystyle x}{\displaystyle l_{(1)(2)}}\right)​=1l(1)(2)[q1⋅(x2−x)+q2⋅(x−x1)]\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right]

Для симплекса s3s_3:

υ(2)(3)(x)\displaystyle \upsilon_{(2)(3)}(x) =q2⋅ϕ2(x)\displaystyle {} = q_2 \cdot \phi_2(x)​+q3⋅ϕ3(x)\displaystyle {} + q_3 \cdot \phi_3(x)​=q2\displaystyle {} = q_2​⋅(x3l(2)(3)−xl(2)(3))\displaystyle {} \cdot \left(\frac{\displaystyle x_3}{\displaystyle l_{(2)(3)}} - \frac{\displaystyle x}{\displaystyle l_{(2)(3)}}\right)​+q3\displaystyle {} + q_3​⋅(−x2l(2)(3)+xl(2)(3))\displaystyle {} \cdot \left(\frac{\displaystyle -x_2}{\displaystyle l_{(2)(3)}} + \frac{\displaystyle x}{\displaystyle l_{(2)(3)}}\right)​=1l(2)(3)[q2⋅(x3−x)+q3⋅(x−x2)]\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right]

Для симплекса s4s_4:

υ(3)(4)(x)\displaystyle \upsilon_{(3)(4)}(x) =q3⋅ϕ3(x)\displaystyle {} = q_3 \cdot \phi_3(x)​+q4⋅ϕ4(x)\displaystyle {} + q_4 \cdot \phi_4(x)​=q3\displaystyle {} = q_3​⋅(x4l(3)(4)−xl(3)(4))\displaystyle {} \cdot \left(\frac{\displaystyle x_4}{\displaystyle l_{(3)(4)}} - \frac{\displaystyle x}{\displaystyle l_{(3)(4)}}\right)​+q4\displaystyle {} + q_4​⋅(−x3l(3)(4)+xl(3)(4))\displaystyle {} \cdot \left(\frac{\displaystyle -x_3}{\displaystyle l_{(3)(4)}} + \frac{\displaystyle x}{\displaystyle l_{(3)(4)}}\right)​=1l(3)(4)[q3⋅(x4−x)+q4⋅(x−x3)]\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} \left[q_3 \cdot (x_4 - x) + q_4 \cdot (x - x_3)\right]

Найдём производные пробных функций по xx:

Для симплекса s1s_1:

υ(0)(1)′(x)\displaystyle \upsilon'_{(0)(1)}(x) =1l(0)(1)(q1−q0)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} (q_1 - q_0)

Для симплекса s2s_2:

υ(1)(2)′(x)\displaystyle \upsilon'_{(1)(2)}(x) =1l(1)(2)(q2−q1)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} (q_2 - q_1)

Для симплекса s3s_3:

υ(2)(3)′(x)\displaystyle \upsilon'_{(2)(3)}(x) =1l(2)(3)(q3−q2)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} (q_3 - q_2)

Для симплекса s4s_4:

υ(3)(4)′(x)\displaystyle \upsilon'_{(3)(4)}(x) =1l(3)(4)(q4−q3)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} (q_4 - q_3)

Вычислим квадраты пробных функций:

Для симплекса s1s_1:

υ(0)(1)2(x)\displaystyle \upsilon_{(0)(1)}^2(x) =1l(0)(1)2[q0⋅(x1−x)+q1⋅(x−x0)]2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} \left[q_0 \cdot (x_1 - x) + q_1 \cdot (x - x_0)\right]^2​=1l(0)(1)2[q02(x1−x)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} \bigl[ q_0^2(x_1 - x)^2​+2q0q1(x1−x)(x−x0)\displaystyle {} + 2q_0q_1(x_1 - x)(x - x_0)​+q12(x−x0)2]\displaystyle {} + q_1^2(x - x_0)^2 \bigr]

Для симплекса s2s_2:

υ(1)(2)2(x)\displaystyle \upsilon_{(1)(2)}^2(x) =1l(1)(2)2[q1⋅(x2−x)+q2⋅(x−x1)]2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right]^2​=1l(1)(2)2[q12(x2−x)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \bigl[ q_1^2(x_2 - x)^2​+2q1q2(x2−x)(x−x1)\displaystyle {} + 2q_1q_2(x_2 - x)(x - x_1)​+q22(x−x1)2]\displaystyle {} + q_2^2(x - x_1)^2 \bigr]

Для симплекса s3s_3:

υ(2)(3)2(x)\displaystyle \upsilon_{(2)(3)}^2(x) =1l(2)(3)2[q2⋅(x3−x)+q3⋅(x−x2)]2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right]^2​=1l(2)(3)2[q22(x3−x)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \bigl[ q_2^2(x_3 - x)^2​+2q2q3(x3−x)(x−x2)\displaystyle {} + 2q_2q_3(x_3 - x)(x - x_2)​+q32(x−x2)2]\displaystyle {} + q_3^2(x - x_2)^2 \bigr]

Для симплекса s4s_4:

υ(3)(4)2(x)\displaystyle \upsilon_{(3)(4)}^2(x) =1l(3)(4)2[q3⋅(x4−x)+q4⋅(x−x3)]2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} \left[q_3 \cdot (x_4 - x) + q_4 \cdot (x - x_3)\right]^2​=1l(3)(4)2[q32(x4−x)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} \bigl[ q_3^2(x_4 - x)^2​+2q3q4(x4−x)(x−x3)\displaystyle {} + 2q_3q_4(x_4 - x)(x - x_3)​+q42(x−x3)2]\displaystyle {} + q_4^2(x - x_3)^2 \bigr]

Вычислим квадраты производных пробных функций:

Для симплекса s1s_1:

(υ(0)(1)′(x))2\displaystyle \left(\upsilon'_{(0)(1)}(x)\right)^2 =1l(0)(1)2(q1−q0)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} (q_1 - q_0)^2​=1l(0)(1)2(q02−2q0q1+q12)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} \left(q_0^2 - 2q_0q_1 + q_1^2\right)

Для симплекса s2s_2:

(υ(1)(2)′(x))2\displaystyle \left(\upsilon'_{(1)(2)}(x)\right)^2 =1l(1)(2)2(q2−q1)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} (q_2 - q_1)^2​=1l(1)(2)2(q12−2q1q2+q22)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left(q_1^2 - 2q_1q_2 + q_2^2\right)

Для симплекса s3s_3:

(υ(2)(3)′(x))2\displaystyle \left(\upsilon'_{(2)(3)}(x)\right)^2 =1l(2)(3)2(q3−q2)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} (q_3 - q_2)^2​=1l(2)(3)2(q22−2q2q3+q32)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \left(q_2^2 - 2q_2q_3 + q_3^2\right)

Для симплекса s4s_4:

(υ(3)(4)′(x))2\displaystyle \left(\upsilon'_{(3)(4)}(x)\right)^2 =1l(3)(4)2(q4−q3)2\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} (q_4 - q_3)^2​=1l(3)(4)2(q32−2q3q4+q42)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} \left(q_3^2 - 2q_3q_4 + q_4^2\right)

Запишем линейную аппроксимацию функции источника f(x)f(x) на каждом симплексе согласно (6.36). В нашем примере источник действует только на внутренних симплексах: на крайних он тождественно равен нулю.

Для симплекса s1s_1 (между x0x_0 и x1x_1):

Поскольку функция плотности тепловых источников тождественно равна нулю на этом интервале:

f~(0)(1)(x)\displaystyle \widetilde{f}_{(0)(1)}(x) =0.\displaystyle {} = 0.

Для симплекса s2s_2 (между x1x_1 и x2x_2):

f~(1)(2)(x)\displaystyle \widetilde{f}_{(1)(2)}(x) =e1(1)(2)\displaystyle {} = e_1^{(1)(2)}​+e2(1)(2)⋅x,\displaystyle {} + e_2^{(1)(2)} \cdot x,

где:

e1(1)(2)\displaystyle e_1^{(1)(2)} =f1⋅x2−f2⋅x1x2−x1\displaystyle {} = \frac{\displaystyle f_1 \cdot x_2 - f_2 \cdot x_1}{\displaystyle x_2 - x_1}​=f1⋅x2−f2⋅x1l(1)(2)\displaystyle {} = \frac{\displaystyle f_1 \cdot x_2 - f_2 \cdot x_1}{\displaystyle l_{(1)(2)}}e2(1)(2)\displaystyle e_2^{(1)(2)} =f2−f1x2−x1\displaystyle {} = \frac{\displaystyle f_2 - f_1}{\displaystyle x_2 - x_1}​=f2−f1l(1)(2)\displaystyle {} = \frac{\displaystyle f_2 - f_1}{\displaystyle l_{(1)(2)}}

Для симплекса s3s_3 (между x2x_2 и x3x_3):

f~(2)(3)(x)\displaystyle \widetilde{f}_{(2)(3)}(x) =e1(2)(3)\displaystyle {} = e_1^{(2)(3)}​+e2(2)(3)⋅x,\displaystyle {} + e_2^{(2)(3)} \cdot x,

где:

e1(2)(3)\displaystyle e_1^{(2)(3)} =f2⋅x3−f3⋅x2x3−x2\displaystyle {} = \frac{\displaystyle f_2 \cdot x_3 - f_3 \cdot x_2}{\displaystyle x_3 - x_2}​=f2⋅x3−f3⋅x2l(2)(3)\displaystyle {} = \frac{\displaystyle f_2 \cdot x_3 - f_3 \cdot x_2}{\displaystyle l_{(2)(3)}}e2(2)(3)\displaystyle e_2^{(2)(3)} =f3−f2x3−x2\displaystyle {} = \frac{\displaystyle f_3 - f_2}{\displaystyle x_3 - x_2}​=f3−f2l(2)(3)\displaystyle {} = \frac{\displaystyle f_3 - f_2}{\displaystyle l_{(2)(3)}}

Для симплекса s4s_4 (между x3x_3 и x4x_4):

Поскольку функция плотности тепловых источников тождественно равна нулю на этом интервале:

f~(3)(4)(x)\displaystyle \widetilde{f}_{(3)(4)}(x) =0.\displaystyle {} = 0.

Теперь запишем произведение функции источника на соответствующие пробные функции для каждого симплекса.

Для симплекса s1s_1:

Поскольку f~(0)(1)(x)\widetilde{f}_{(0)(1)}(x)​=0{} = 0, имеем:

f~(0)(1)(x)⋅υ(0)(1)(x)\displaystyle \widetilde{f}_{(0)(1)}(x) \cdot \upsilon_{(0)(1)}(x) =0.\displaystyle {} = 0.

Для симплекса s2s_2:

f~(1)(2)(x)⋅υ(1)(2)(x)\displaystyle \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) =(e1(1)(2)+e2(1)(2)⋅x)\displaystyle {} = \left(e_1^{(1)(2)} + e_2^{(1)(2)} \cdot x\right) \quad​⋅1l(1)(2)[q1⋅(x2−x)+q2⋅(x−x1)].\displaystyle {} \cdot \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right].

Подставляя значения e1(1)(2)e_1^{(1)(2)} и e2(1)(2)e_2^{(1)(2)}:

f~(1)(2)(x)⋅υ(1)(2)(x)\displaystyle \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) =1l(1)(2)2[(f1⋅x2−f2⋅x1)+(f2−f1)⋅x]\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left[(f_1 \cdot x_2 - f_2 \cdot x_1) + (f_2 - f_1) \cdot x\right] \quad​⋅[q1⋅(x2−x)+q2⋅(x−x1)].\displaystyle {} \cdot \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right].

Преобразуем первый множитель:

(f1⋅x2−f2⋅x1)\displaystyle (f_1 \cdot x_2 - f_2 \cdot x_1)​+(f2−f1)⋅x\displaystyle {} + (f_2 - f_1) \cdot x​=f1⋅(x2−x)\displaystyle {} = f_1 \cdot (x_2 - x)​+f2⋅(x−x1).\displaystyle {} + f_2 \cdot (x - x_1).

Следовательно:

f~(1)(2)(x)⋅υ(1)(2)(x)\displaystyle \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) =1l(1)(2)2[f1⋅(x2−x)+f2⋅(x−x1)]\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left[f_1 \cdot (x_2 - x) + f_2 \cdot (x - x_1)\right] \quad​⋅[q1⋅(x2−x)+q2⋅(x−x1)].\displaystyle {} \cdot \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right].

Раскроем скобки:

f~(1)(2)(x)⋅υ(1)(2)(x)\displaystyle \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) =1l(1)(2)2[f1q1(x2−x)2+f1q2(x2−x)(x−x1)+f2q1(x−x1)(x2−x)+f2q2(x−x1)2].\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \Big[ f_1 q_1 (x_2 - x)^2 + f_1 q_2 (x_2 - x)(x - x_1) + f_2 q_1 (x - x_1)(x_2 - x) + f_2 q_2 (x - x_1)^2 \Big].

Для симплекса s3s_3:

f~(2)(3)(x)⋅υ(2)(3)(x)\displaystyle \widetilde{f}_{(2)(3)}(x) \cdot \upsilon_{(2)(3)}(x) =(e1(2)(3)+e2(2)(3)⋅x)\displaystyle {} = \left(e_1^{(2)(3)} + e_2^{(2)(3)} \cdot x\right) \quad​⋅1l(2)(3)[q2⋅(x3−x)+q3⋅(x−x2)].\displaystyle {} \cdot \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right].

По аналогии с симплексом s2s_2:

f~(2)(3)(x)⋅υ(2)(3)(x)\displaystyle \widetilde{f}_{(2)(3)}(x) \cdot \upsilon_{(2)(3)}(x) =1l(2)(3)2[f2⋅(x3−x)+f3⋅(x−x2)]\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \left[f_2 \cdot (x_3 - x) + f_3 \cdot (x - x_2)\right] \quad​⋅[q2⋅(x3−x)+q3⋅(x−x2)].\displaystyle {} \cdot \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right].

Раскроем скобки:

f~(2)(3)(x)⋅υ(2)(3)(x)\displaystyle \widetilde{f}_{(2)(3)}(x) \cdot \upsilon_{(2)(3)}(x) =1l(2)(3)2[f2q2(x3−x)2+f2q3(x3−x)(x−x2)+f3q2(x−x2)(x3−x)+f3q3(x−x2)2].\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \Big[ f_2 q_2 (x_3 - x)^2 + f_2 q_3 (x_3 - x)(x - x_2) + f_3 q_2 (x - x_2)(x_3 - x) + f_3 q_3 (x - x_2)^2 \Big].

Для симплекса s4s_4:

Поскольку f~(3)(4)(x)\widetilde{f}_{(3)(4)}(x)​=0{} = 0, имеем:

f~(3)(4)(x)⋅υ(3)(4)(x)\displaystyle \widetilde{f}_{(3)(4)}(x) \cdot \upsilon_{(3)(4)}(x) =0.\displaystyle {} = 0.

Перейдём к граничному члену. Граница одномерной области — две точки x0x_0 и x4x_4, поэтому поверхностный интеграл сводится к сумме по граничным узлам:

∫Sυ⋅∇υ dS\displaystyle \int_S \upsilon \cdot \nabla \upsilon \,dS =q0⋅∂u∂n∣x=x0\displaystyle {} = q_0 \cdot \frac{\partial u}{\partial n}\bigg|_{x=x_0}​+q4⋅∂u∂n∣x=x4.\displaystyle {} + q_4 \cdot \frac{\partial u}{\partial n}\bigg|_{x=x_4}.

Поскольку в этом примере заданы условия Неймана, нормальная производная на границе известна: ∂u∂n∣x=x0\frac{\partial u}{\partial n}\bigg|_{x=x_0}​=g0{} = g_0, ∂u∂n∣x=x4\frac{\partial u}{\partial n}\bigg|_{x=x_4}​=g4{} = g_4. Этот член входит только в уравнения граничных узлов и переносится в правую часть как вектор потока; в матрицу системы он не входит.

Интеграл жёсткости для каждого симплекса вычисляется как интеграл квадрата производной пробной функции. Поскольку производные пробных функций являются константами, интегрирование сводится к умножению на длину симплекса.

Согласно локальной формуле из раздела «Матрица жёсткости 1D», для каждого симплекса этот интеграл равен 1l(i)(i+1)(qi2−2qiqi+1+qi+12)\frac{1}{l_{(i)(i+1)}}\left(q_i^2 - 2 q_i q_{i+1} + q_{i+1}^2\right).

Общий интеграл жёсткости по всей области MM равен сумме интегралов по всем симплексам:

∫M(∇υ)2 dM\displaystyle \int_M (\nabla \upsilon)^2 \,dM =∑i=14∫si(υ′(x))2 dx\displaystyle {} = \sum_{i=1}^{4} \int_{s_i} \left(\upsilon'(x)\right)^2 \,dx=1l(0)(1)(q02−2q0q1+q12)\displaystyle {} = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} \left(q_0^2 - 2q_0q_1 + q_1^2\right)​+1l(1)(2)(q12−2q1q2+q22)\displaystyle {} + \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} \left(q_1^2 - 2q_1q_2 + q_2^2\right)\displaystyle \quad​+1l(2)(3)(q22−2q2q3+q32)\displaystyle {} + \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} \left(q_2^2 - 2q_2q_3 + q_3^2\right)​+1l(3)(4)(q32−2q3q4+q42).\displaystyle {} + \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} \left(q_3^2 - 2q_3q_4 + q_4^2\right).

Интеграл демпфирования для каждого симплекса вычисляется как интеграл квадрата пробной функции.

Согласно локальной формуле из раздела «Матрица демпфирования 1D», интеграл квадрата пробной функции на симплексе равен l(i)(i+1)3(qi2+qiqi+1+qi+12)\frac{l_{(i)(i+1)}}{3}\left(q_i^2 + q_i q_{i+1} + q_{i+1}^2\right).

Общий интеграл демпфирования по всей области MM равен сумме интегралов по всем симплексам:

∫Mυ2 dM\displaystyle \int_M \upsilon^2 \,dM =∑i=14∫siυ2(x) dx\displaystyle {} = \sum_{i=1}^{4} \int_{s_i} \upsilon^2(x) \,dx=l(0)(1)3(q02+q0q1+q12)\displaystyle {} = \frac{\displaystyle l_{(0)(1)}}{\displaystyle 3} \left(q_0^2 + q_0q_1 + q_1^2\right)​+l(1)(2)3(q12+q1q2+q22)\displaystyle {} + \frac{\displaystyle l_{(1)(2)}}{\displaystyle 3} \left(q_1^2 + q_1q_2 + q_2^2\right)\displaystyle \quad​+l(2)(3)3(q22+q2q3+q32)\displaystyle {} + \frac{\displaystyle l_{(2)(3)}}{\displaystyle 3} \left(q_2^2 + q_2q_3 + q_3^2\right)​+l(3)(4)3(q32+q3q4+q42).\displaystyle {} + \frac{\displaystyle l_{(3)(4)}}{\displaystyle 3} \left(q_3^2 + q_3q_4 + q_4^2\right).

Интеграл нагрузки для каждого симплекса вычисляется как интеграл произведения функции источника на пробную функцию.

Согласно локальной формуле из раздела «Вектор нагрузки 1D», на крайних симплексах s1s_1 и s4s_4 интеграл равен нулю (источник там тождественно нулевой), а на s2s_2 и s3s_3 — соответствующему вкладу локального вектора нагрузки.

Общий интеграл нагрузки по всей области MM равен сумме интегралов по всем симплексам:

∫Mf~(x)⋅υ(x) dM\displaystyle \int_M \widetilde{f}(x) \cdot \upsilon(x) \,dM =∑i=14∫sif~(x)⋅υ(x) dx\displaystyle {} = \sum_{i=1}^{4} \int_{s_i} \widetilde{f}(x) \cdot \upsilon(x) \,dx=0\displaystyle {} = 0​+l(1)(2)6(2f1q1+f1q2+f2q1+2f2q2)\displaystyle {} + \frac{\displaystyle l_{(1)(2)}}{\displaystyle 6} \left(2f_1 q_1 + f_1 q_2 + f_2 q_1 + 2f_2 q_2\right)\displaystyle \quad​+l(2)(3)6(2f2q2+f2q3+f3q2+2f3q3)\displaystyle {} + \frac{\displaystyle l_{(2)(3)}}{\displaystyle 6} \left(2f_2 q_2 + f_2 q_3 + f_3 q_2 + 2f_3 q_3\right)​+0.\displaystyle {} + 0.

Упростив:

∫Mf~(x)⋅υ(x) dM\displaystyle \int_M \widetilde{f}(x) \cdot \upsilon(x) \,dM =l(1)(2)6(2f1q1+f1q2+f2q1+2f2q2)\displaystyle {} = \frac{\displaystyle l_{(1)(2)}}{\displaystyle 6} \left(2f_1 q_1 + f_1 q_2 + f_2 q_1 + 2f_2 q_2\right) \quad​+l(2)(3)6(2f2q2+f2q3+f3q2+2f3q3).\displaystyle {} + \frac{\displaystyle l_{(2)(3)}}{\displaystyle 6} \left(2f_2 q_2 + f_2 q_3 + f_3 q_2 + 2f_3 q_3\right).

Запишем все полученные интегралы в матричной форме. Введём вектор неизвестных q→\overrightarrow{q}​=(q0q1q2q3q4)T{} = \begin{pmatrix} q_0 & q_1 & q_2 & q_3 & q_4 \end{pmatrix}^T.

Интеграл жёсткости можно записать в матричном виде:

∫M(∇υ)2 dM\displaystyle \int_M (\nabla \upsilon)^2 \,dM =q→T⋅K⋅q→,\displaystyle {} = \overrightarrow{q}^T \cdot K \cdot \overrightarrow{q},

где KK — матрица жёсткости размером 5×55 \times 5:

K\displaystyle K =(1l(0)(1)−1l(0)(1)000−1l(0)(1)1l(0)(1)+1l(1)(2)−1l(1)(2)000−1l(1)(2)1l(1)(2)+1l(2)(3)−1l(2)(3)000−1l(2)(3)1l(2)(3)+1l(3)(4)−1l(3)(4)000−1l(3)(4)1l(3)(4)).\displaystyle {} = \begin{pmatrix} \frac{1}{l_{(0)(1)}} & -\frac{1}{l_{(0)(1)}} & 0 & 0 & 0 \\ -\frac{1}{l_{(0)(1)}} & \frac{1}{l_{(0)(1)}} + \frac{1}{l_{(1)(2)}} & -\frac{1}{l_{(1)(2)}} & 0 & 0 \\ 0 & -\frac{1}{l_{(1)(2)}} & \frac{1}{l_{(1)(2)}} + \frac{1}{l_{(2)(3)}} & -\frac{1}{l_{(2)(3)}} & 0 \\ 0 & 0 & -\frac{1}{l_{(2)(3)}} & \frac{1}{l_{(2)(3)}} + \frac{1}{l_{(3)(4)}} & -\frac{1}{l_{(3)(4)}} \\ 0 & 0 & 0 & -\frac{1}{l_{(3)(4)}} & \frac{1}{l_{(3)(4)}} \end{pmatrix}.
(J.1)

Матрица KK является симметричной и трёхдиагональной.

Интеграл демпфирования можно записать в матричном виде:

∫Mυ2 dM\displaystyle \int_M \upsilon^2 \,dM =q→T⋅D⋅q→,\displaystyle {} = \overrightarrow{q}^T \cdot D \cdot \overrightarrow{q},

где DD — матрица демпфирования размером 5×55 \times 5:

D\displaystyle D =(l(0)(1)3l(0)(1)6000l(0)(1)6l(0)(1)3+l(1)(2)3l(1)(2)6000l(1)(2)6l(1)(2)3+l(2)(3)3l(2)(3)6000l(2)(3)6l(2)(3)3+l(3)(4)3l(3)(4)6000l(3)(4)6l(3)(4)3).\displaystyle {} = \begin{pmatrix} \frac{l_{(0)(1)}}{3} & \frac{l_{(0)(1)}}{6} & 0 & 0 & 0 \\ \frac{l_{(0)(1)}}{6} & \frac{l_{(0)(1)}}{3} + \frac{l_{(1)(2)}}{3} & \frac{l_{(1)(2)}}{6} & 0 & 0 \\ 0 & \frac{l_{(1)(2)}}{6} & \frac{l_{(1)(2)}}{3} + \frac{l_{(2)(3)}}{3} & \frac{l_{(2)(3)}}{6} & 0 \\ 0 & 0 & \frac{l_{(2)(3)}}{6} & \frac{l_{(2)(3)}}{3} + \frac{l_{(3)(4)}}{3} & \frac{l_{(3)(4)}}{6} \\ 0 & 0 & 0 & \frac{l_{(3)(4)}}{6} & \frac{l_{(3)(4)}}{3} \end{pmatrix}.

Матрица DD также является симметричной и трёхдиагональной.

Интеграл нагрузки можно записать в матричном виде:

∫Mf~(x)⋅υ(x) dM\displaystyle \int_M \widetilde{f}(x) \cdot \upsilon(x) \,dM =F→T⋅q→,\displaystyle {} = \overrightarrow{F}^T \cdot \overrightarrow{q},

где F→\overrightarrow{F} — вектор нагрузки размером 5×15 \times 1:

F→\displaystyle \overrightarrow{F} =(02f1⋅l(1)(2)6+f2⋅l(1)(2)6f1⋅l(1)(2)6+2f2⋅l(1)(2)6+2f2⋅l(2)(3)6+f3⋅l(2)(3)6f2⋅l(2)(3)6+2f3⋅l(2)(3)60).\displaystyle {} = \begin{pmatrix} 0 \\ \frac{2f_1 \cdot l_{(1)(2)}}{6} + \frac{f_2 \cdot l_{(1)(2)}}{6} \\ \frac{f_1 \cdot l_{(1)(2)}}{6} + \frac{2f_2 \cdot l_{(1)(2)}}{6} + \frac{2f_2 \cdot l_{(2)(3)}}{6} + \frac{f_3 \cdot l_{(2)(3)}}{6} \\ \frac{f_2 \cdot l_{(2)(3)}}{6} + \frac{2f_3 \cdot l_{(2)(3)}}{6} \\ 0 \end{pmatrix}.

Упростив, получим:

F→\displaystyle \overrightarrow{F} =(0l(1)(2)6(2f1+f2)f1⋅l(1)(2)6+f2(l(1)(2)+l(2)(3))3+f3⋅l(2)(3)6l(2)(3)6(f2+2f3)0).\displaystyle {} = \begin{pmatrix} 0 \\ \frac{l_{(1)(2)}}{6} \left(2f_1 + f_2\right) \\ \frac{f_1 \cdot l_{(1)(2)}}{6} + \frac{f_2 (l_{(1)(2)} + l_{(2)(3)})}{3} + \frac{f_3 \cdot l_{(2)(3)}}{6} \\ \frac{l_{(2)(3)}}{6} \left(f_2 + 2f_3\right) \\ 0 \end{pmatrix}.

Граничный интеграл, как показано выше, в матрицу не сворачивается — он даёт вектор потока в правую часть, ненулевой только в граничных узлах:

a2s→\displaystyle a^2 \overrightarrow{s} =a2\displaystyle {} = a^2​(g0000g4).\displaystyle \begin{pmatrix} g_0 \\ 0 \\ 0 \\ 0 \\ g_4 \end{pmatrix}.

Стационарное уравнение

Согласно (I.8), стационарная система имеет вид

a2K⋅q→\displaystyle a^2 K \cdot \overrightarrow{q} =F→\displaystyle {} = \overrightarrow{F}​+a2s→.\displaystyle {} + a^2 \overrightarrow{s}.

Нестационарное уравнение

Согласно (I.9), нестационарная система:

[D+Δt⋅a2K]⋅q→n\displaystyle \left[D + \Delta t \cdot a^2 K\right] \cdot \overrightarrow{q}_n =Δt⋅F→\displaystyle {} = \Delta t \cdot \overrightarrow{F}​+D⋅q→n−1\displaystyle {} + D \cdot \overrightarrow{q}_{n-1}​+Δt⋅a2s→.\displaystyle {} + \Delta t \cdot a^2 \overrightarrow{s}.

Обозначения:

  • KK — матрица жёсткости 5×55 \times 5,
  • DD — матрица демпфирования 5×55 \times 5,
  • F→\overrightarrow{F} — вектор нагрузки 5×15 \times 1,
  • q→\overrightarrow{q} — вектор неизвестных 5×15 \times 1 для стационарного случая,
  • q→n,q→n−1\overrightarrow{q}_n, \overrightarrow{q}_{n-1} — векторы неизвестных на текущем и предыдущем временных слоях,
  • aa — коэффициент температуропроводности,
  • Δt\Delta t — временной шаг.

Матрицы KK и DD симметричны и разрежены (трёхдиагональны). Оператор стационарной системы — полная матрица жёсткости a2Ka^2 K; для чисто неймановой задачи она вырождена (K⋅1→K\cdot\overrightarrow{1}​=0{} =0), поэтому решение определено с точностью до аддитивной константы и замыкается условием нормировки. Подробнее учёт условий Неймана разобран в приложении «Учёт граничных условий Неймана».