Метод сеток для уравнений Пуассона и Лапласа

 

Рассмотрим первую краевую задачу для уравнения Пуассона в прямоугольной области 0 ≤ x ≤ α; 0 ≤ y ≤ b;

U_xx+U_yy=f (x,y)  ,                                        (12.7)

 (ГУ) U_r= φ (x,y)  .                                       (12.8)

 

Выбрав шаги по времени и по координатам, строим сетку

x_l = l · h; x ₀ = 0; x_l = l · h; x_n = α; y_J = J · τ; y ₀ = 0; y_m = b

где         t = 1, 2 ,….,n; J = 1, 2 ,……., m.

Заменяя в каждом внутреннем узле (x_l, y_J) производные конечными разностями по формулам (12.4), получаем конечно-разностные уравнения

 ,                 (12.9)

U _{0, J} = φ (0, y_J) = φ _{0, J}; U_{n , J} = φ (α, y_J) = φ_{n , J},                    (12.10)

U_{l , 0} = φ (x, 0) = φ_{l , 0}; U_{l , m} = φ (x, b) = φ_{l , m},                         (12.11)

 

где t = 1, 2 ,…., n; J = 1, 2 ,……., m.

Уравнения (12.9) вместе с условиями (12.10)-(12.11) образуют систему линейных алгебраических уравнений относительно U_{l , J}. Число неизвестных и число уравнений в системе равно числу внутренних узлов сетки. Схема узловых точек для уравнений (12.9) изображена на рис. 12.1.

Если   h = τ    и J_{l , J} = 0, то конечно-разностные уравнения имеют вид

 

             (12.12)

 

Пример 1. Найти решение уравнения теплопроводности

U_xx = U_t; 0 ≤ x ≤ 1; t ≥ 0;

                        (НУ) U (x, 0) =4 · x · (1 -x);

(ГУ) U (0 ,t) =U (1 ,t) = 0

методом сеток с точностью до ε = 0.001 на первых трех временных уровнях, если шаг по координате Δ x = h = 0,1, по времени Δ t = τ = 1/600.

 

 

Решение. Для вычисления U_{l , J} воспользуемся системой конечно-разностных уравнений (12.9) при g = 0:

=  ( + )+  ,              (12.13)

(НУ) U_{l , 0} = 4 · x_l · (1- x_l); l = 1, 2 ,…., 10,                     (12.14)

 (ГУ) U _{0, J} = U _{10, J} = 0; J = 0,1,2.                                 (12.15)

 

Составим таблицу значений U_{l , J}  (табл. 12.1).

 

Таблица 12.1

l J	0	1	2	3	4	5	6	7	8	9	10
0	0	0.360	0.640	0.840	0.960	1	0.960	0.840	0.640	0.360	0
1	0	0.347	0.627	0.827	0.977	0.987	0.977	0.827	0.627	0.347	0
2	0	0.336	0.613	0.813	0.933	0.973	0.933	0.813	0.613	0.336	0

 

В силу симметрии начальных  и краевых условий относительно оси   x = 0,5 решение U_{l , J} будет симметричным относительно этой же оси.

Это позволяет вести расчеты для t = 1, 2 ,…., 5, а для остальных значений следует заполнить таблицу, исходя из свойства симметрии:

 

U _{0, J} = U _{10, J}; U _{1, J} = U _{9, J}; U _{2, J} = U _{8, J}; U _{3, J} = U _{7, J}; U _{4, J} = U _{6, J}   .

Начальная строка этой таблицы (J = 0) заполняется на основании зaданного начального условия (12.14).

В первый (l = 0) и последний (l = 10) столбцы вписываются данные граничных условий (12.15). Остальные строки таблицы последовательно заполняются с помощью применения расчетной формулы (12.13). На первом шаге

 

U_{l , J} = (1/6) ·  .

Например,

U _1,1 = (1/6) ·  = (1/6) · (0,640+4∙0,360) = 0,347.

На втором шаге

U_{l ,2} = (1/6) ·  .

 

 

Пример 2. Применяя метод сеток с шагом h = τ = 1/3, найти решениеуравнения Лапласа в квадрате с вершинами А (0,0), В (0,1), С (1,1), D (1,0), удовлетворяющее краевым условиям

U (0 ,y) = 30 · y; U (x, 1) = 30 · (1 -x ²); U (1 ,y) = 0; U (x, 0) = 0.

Решение. Конечно-разностные уравнения для уравнения Лапласа при                h = τ  = 1/3имеют вид:

U _{0, J} = 0,25 · ,               (12.16)

     (ГУ)    U _{0, J} = 30 · J · h; U _{3, J} = 0 ;                                   (12.17)

 (НУ)  (U_{l, 0} = 0; U_{l, 3} = 30 · (1 - (l · h)²),                               (12.18)

где x = l · h, y = J · h, l = 1,2,3, J = 0,1,2,3.

Составим таблицу U_{l , J} (табл. 12.2)

Таблица 12.2

y x	0	1/3	2/3	1
0	U 0,0 = 0	U0 ,1 = 10	U 0,2 = 20	U 0,3 = 30
1/3	U l,0 = 0	U l,1	U l,2	U l,3 = 26,67
2/3	U 2,0 = 0	U 2,1	U 2,2	U 2,3 = 16,67
1	U 3,0 = 0	U 3,1 = 0	U 3,2 = 0	U 3,3 = 0

 

Первая и последняя строки этой таблицы получены из граничных условий (12.17), а первый и последний столбцы - из (12.18).

 

U _l,1, U _l,2, U _2,1, U _2,2 находимиз следующей системы линейных алгебраических уравнений, которые получены из расчетной формул (12.16):

 

 

.

 

 

ОБРАЗЦЫ ВЫПОЛНЕНИЯ КОНТРОЛЬНЫХ РАБОТ

Контрольная работа № 5

Вычислить значение интеграла

, , а = 0, b = 0, .

Выберем формулу для приближенного вычисления заданного определенного интеграла. Для чего найдем по формулам (9.4), (9.6), (9.10) число n точек разбиения отрезка [0,1] на частичные, которые обеспечат требуемую точность при вычислении по формулам прямоугольников, трапеций и парабол соответственно. А затем остановимся на той из приближенных формул, для которой число n будет наименьшим.

Чтобы воспользоваться формулами (9.4), (9.6), (9.10), вычислим и оценим первую, вторую и четвертую производные подынтегральной функции  на отрезке [0,1]:

, ,

,

и так как функция f (x) и ее производные убывают на отрезке [0,1], то

,

,

.

Найдем n. Для формулы прямоугольников из (9.4) получаем

.

Для формулы трапеций из (9.6) получаем

, .

Для формулы парабол из (9.10) получаем

; .

 

Таким образом, наименьшего объема вычислений при одинаковой точности потребует формула парабол (9.8) – n = 2 m = 8 (n должно быть четным), применяя которую и вычислим приближенно заданный интеграл.

По числу n = 2 m = 8 найдем шаг интегрирования .

Составим таблицу (табл.1.) значений подынтегральной функции в точках x_i = ih, , записывая ординаты с четными и нечетными номерами в разные столбцы. В последней строке таблицы запишем результаты суммирования по этим столбцам. Вычисление будем вести с четырьмя знаками после запятой, а окончательный ответ округлим до трех знаков после запятой.

Применяя формулу парабол (9.8), получим

Вычислим заданный интеграл по формуле Ньютона–Лейбница

.

Итак, требуемая точность вычислений достигнута.

 

Таблица 1

i	xi = a+ih,	1+2х i	Значения у i =1/(1+2х i)
			При i =0, i =8	При четном i	При нечетном i
0	0	1,0
1.	0,125	1,250			0,80
2.	0,25	1,50		0,6667
3.	0,375	1,750			0,5714
4.	0,50	2,00		0,50
5.	0,625	2,250			0,4444
6.	0,750	2,50		0,40
7.	0,875	2,750			0,3636
8.	1,0	3,0	0,3333

Суммы