Решение дифференциальных  уравнений с частными  производными

Пусть требуется решить смешанную задачу для уравнения теплопроводности, т.е. найти функцию U (x, t), удовлетворяющую дифференциальному уравнению с частными производными

_,                                                                     (18)

где U = U (x, t) – температура в сечении стержня с абсциссой х в момент времени t, если известны:

температура в сечении стержня с абсциссой х в момент времени t =0:

U (х, 0) = f (x) для 0 £ x £ L,                          (19)

и температура на концах стержня в любой момент времени t:

U (0, t) = α (t), U (L, t) = β (t) для 0 £ t £ T,            (20)

где f (x), α (t), β (t) – известные непрерывные функции.

Будем предполагать, что смешанная задача (18)-(20) имеет единственное решение U = U (x, t) в области  

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

Введем прямоугольную сетку, построенную на области D:

 x_i = ih для i = 0, 1, 2,….,n,

 t_k = kd для k = 0,1, 2,….,m,

где h  =  – шаг сетки по х,  d =  – шаг сетки по t. При этом прямоугольная область D будет разбита на прямоугольные части прямыми х = x_i  и t = t_k  в системе координат Охt. Точки (x_i, t_k),лежащие на пересечении этих прямых, называют узлами двумерной сетки, или просто узлами.

Вместо неизвестной функции U (x, t) будем рассматривать функцию, заданную в узлах сетки: u_{i, k} = u(i, k), которую называют сеточной функцией. Будем искать значения сеточной функции, соответствующие приближенному решению задачи, т.е. u_{i, k} » U (x_i, t_k).

Неизвестные u_{0, k}  и u_{n, k}  для k = 0, 1,….,m находят из граничных условий (20):

u_{0, k} = U (0, t_k) = α (t_k), u_{n, k} = U (L, t_k) = β (t_k).

Значения u_i,0, соответствующие моменту времени t₀ = 0, находят из начальных условий (19):

u_{i, 0} = U (x_i, 0) = f (x_i), для i= 1,….,n –1.

Для определения значений u_{i, k} во внутренних узлах заменим частные производные отношениями конечных разностей по приближенным формулам:

,

тогда уравнению (18) будет соответствовать сеточное уравнение:

.

Если выбрать , то получим простую формулу для вычисления значений сеточной функции: .

     Можно доказать, что для значений  получается устойчивое решение смешанной задачи, т.е. малым изменениям исходных данных задачи будут соответствовать лишь малые изменения полученного решения.

Если выбрать шаги сетки h и d таким образом, чтобы было выполнено условие , то получается наиболее точная формула для расчетов:

,

для i = 1, 2,….,n – 1, k = 0, 1,….,m – 1.

Значения u_{i, k} для i = 1, 2,….,n – 1 находят «по слоям»: сначала вычисляют u_{i 0} на «нулевом слое», затем вычисляют u_i,1, u_i,2 ,… – значения на каждом следующем, «(k +1)-м слое», соответствующем моменту времени t_k+1, используя значения на предыдущем, «k-м слое».

Окончательно получаем разностную модель задачи – расчетные формулы для вычисления всех значений функции U (x_i, t_k) = u_{i, k} в узлах сетки, покрывающей область D:

,                               (21)

,                               (22)

,                             (23)

(24)

Решение примерного варианта контрольной работы

Задача 1. Даны  приближенные  значения  величин  x = 1,6,  y = 0,35,   где D(х)= 0,02, d (y)= 4%. Требуется:

1) вычислить значение величины , оценить предельную абсолютную погрешность D(s) и округлить значение s в соответствии с погрешностью;

2) вычислить приближенное значение функции f (x, y) = ·(2x – 3y), оценить предельную абсолютную погрешность значения функции и округлить его в соответствии с погрешностью.

Решение.

1). Вычислим приближенное значение s: s = .

Оценим предельную абсолютную погрешность результата по формуле (2), в которой

D(х)= 0,02; D(y) = y·d (y) = 0,35·0,04 = 0,014,

следовательно, по формуле (1) получаем:

D(3y – х) = 3D(y) + D(х) = 0,062.

Предельная абсолютная погрешность значения s:

.

Округлим значение s, оставляя столько же цифр после запятой, сколько их в записи абсолютной погрешности результата: D(s) = 0,1, следовательно, значение s округляем до одного знака после запятой: s » – 0,6.

2). Вычислим приближенное значение функции f (x, y) = ·(2x – 3y):

f (1,6; 0,35) = ·(3,2 – 1,05) » 1,2649·2,15 » 2,7196.

Оценим предельную абсолютную погрешность результата по формуле (4), в которой:

;

Максимальные значения модулей частных производных:

.

Предельная абсолютная погрешность значения функции:

(округление абсолютной погрешности производится «с избытком»).

Округлим значение функции, оставляя столько же цифр после запятой, сколько их в записи абсолютной погрешности результата: D(f) = 0,2, тогда f (1,6; 0,35) » 2,7.

Ответы: 1) s » – 0,6; D(s) = 0,1;

2) f (1,6; 0,35) » 2,7; D(f) = 0,2.

Задача 2. Дано уравнение: x⁵ + 3x – 2,5 = 0. Требуется:

1)  определить число корней уравнения и найти промежутки их изоляции;

2)  вычислить значение одного из корней уравнения с точностью ε = 0,01 при помощи метода деления отрезка пополам.

Решение.

1). Будем искать такой интервал (a, b), в котором содержится корень уравнения x⁵ + 3x – 2,5 = 0, причем только один. Для этого найдем производную  = 5х⁴ + 3. Из того, что > 0 для любых значений x, следует, что функция f (x) = x⁵ + 3x – 2,5 монотонно возрастающая, значит, уравнение f (x) = 0 имеет только один корень (точку пересечения графика y = f (x) с осью абсцисс).

Подбором находим две точки, в которых функция имеет разные знаки: например, точки a = 0,5, b = 1 (чем ближе друг к другу эти точки, тем меньше потребуется сделать вычислений на этапе уточнения корня). Проверка: , следовательно,  и корень уравнения .

     2). Уточним корень уравнения x⁵ + 3x – 2,5 = 0 с точностью ε = 0,01 методом деления отрезка пополам.

Для этого обозначим: номер шага k = 0; a₀ = a = 0,5; b₀= b = 1, и вычислим начальное значение корня =0,75 и значение функции f (x₀) = + 3x₀–2,5 » – 0,0127 (все промежуточные вычисления будем производить, используя 4 десятичных знака после запятой).

Для k = 1 находим a_k, b_k. По формулам (8.1) –(8.2). Поскольку f (x₀) < 0, то выполнены условия , , поэтому a₁:= x₀ и b₁:= b₀, т.е. [a₁;b₁] = [0,75; 1]. После этого вычислим =0,875, значение функции f (x₁) »0,6379 и проверим выполнение условия (9): . Получаем: . Условие не выполнено, следовательно, переходим к следующему шагу.

Для k = 2, 3,…. находим a_k, b_k, x_k по формулам (8.1) – (8.3), проверяя на каждом шагу выполнение условия (9):  и, если оно выполнено, определяем знак f (x_k). Если f (x_k) < 0, то на следующем шаге изменяем a_k, а если f (x_k) > 0, то изменяем b_k, чтобы на каждом шаге выполнялось условие . Результаты расчетов заносим в таблицу:

k	ak	bk			f (xk) = + 3xk–2,5
0	0,5	1	0,75	0,5	– 0,0127 < 0
1	0,75	1	0,875	0,25	0,6379 > 0
2	0,75	0,875	0,8125	0,125	0,2916 > 0
3	0,75	0,8125	0,7813	0,0625	0,1348 > 0
4	0,75	0,7813	0,7656	0,0313	0,0600 > 0
5	0,75	0,7656	0,7578	0,0156

Процесс вычислений закончен, т.к. . Последнее значение x₅ округляем в соответствии с заданной точностью ε = 0,01:

х* » x₅ » 0,76.

Ответ: .

Задача 3. Дана система линейных алгебраических уравнений:

Требуется:

1) решить систему методом Гаусса;

2) вычислить определитель матрицы системы, используя метод Гаусса.

Решение.

1) Решение системы.

Прямой ход. Построим расширенную матрицу (A|В) размерности n (n + 1), т.е. 5 6, присоединив к матрице А столбец свободных членов В:

Приведение системы к верхнему треугольному виду осуществляем за 3 шага, используя формулы (10), (11).

1шаг (k = 1): исключаем неизвестную х₁ из уравнений с номерами i =2, 3, 4.

Все элементы 1-й (ведущей) строки расширенной матрицы делим на число d₁ = а₁₁ = 2, получаем строку:  = (1; 1,5; –0,5; –2,5; –0,5).

Из 2-й строки вычитаем строку , умноженную на число с₂₁ = а₂₁= 1:

= (1; –2; 3; 4; 1) – (1; 1,5; –0,5; –2,5; –0,5) = (0; –3,5; 3,5; 6,5; 1,5).

Из 3-й строки вычитаем строку , умноженную на число с₃₁ = а₃₁= 5, получаем:

= (5; 1; –1; –8; –1) – 5 (1; 1,5; –0,5; –2,5; –0,5) = (0; –6,5; 1,5; 4,5; 1,5).

Из 4-й строки вычитаем строку , умноженную на число с₄₁ = а₄₁= 1:

= (1; –5; 2; 3; 7) – (1; 1,5; –0,5; –2,5; –0,5) = (0; –6,5; 2,5; 5,5; 7,5).

В результате 1-го шага произошло исключение неизвестной х₁ из 2-го, 3-го и 4-го уравнений системы, расширенная матрица и система теперь имеют вид:

(A|В)′ = ,

    Далее работаем только с 2, 3, 4 строками расширенной матрицы (A|B)′.

2шаг (k = 2): исключаем неизвестную х₂ из уравнений с номерами i =3,4.

Все элементы 2-й (ведущей) строки делим на число d₂ = a′₂₂ = –3,5, получаем строку:  » (0; 1; –1; –1,8571; –0,4286) (все промежуточные вычисления будем производить, округляя результаты до 4 десятичных знаков после запятой).

Из 3-й строки вычитаем строку , умноженную на с₃₂ = a′₃₂ = –6,5:

 » (0; –6,5; 1,5; 4,5; 1,5) + 6,5 (0; 1; –1; –1,8571; –0,4286) »

» (0; 0; –5; –7,5712; –1,2859).

Из 4-й строки вычитаем строку , умноженную на с₄₂ = a′₄₂ = –6,5:

 » (0; –6,5; 2,5; 5,5; 7,5) + 6,5 (0; 1; –1; –1,8571; –0,4286) »

» (0; 0; –4; –6,5712; 4,7141).

В результате 2-го шага произошло исключение неизвестной х₂ из 3-го и 4-го уравнений системы, расширенная матрица и система теперь имеют вид:

(A|В)′′= ,

Далее работаем только с 3, 4 строками расширенной матрицы (A|B)′′.

3шаг (k = 3): исключаем неизвестную х₃ из уравнения с номером i =4.

Все элементы 3-й (ведущей) строки делим на число d₃ = a′′₃₃ = –5, получаем строку: » (0; 0; 1; 1,5142; 0,2572).

Из 4-й строки вычитаем строку , умноженную на с₄₃ = a′′₄₃ = – 4:

» (0; 0; –4; –6,5712; 4,7141) + 4 (0; 0; 1;1,5142; 0,2572)»

» (0; 0; 0; –0,5144; 5,7429).

В результате 3-го шага произошло исключение неизвестной х₃ из 4-го уравнения системы, расширенная матрица теперь имеет вид:

(A|В)′′′= ,

Таким образом, в результате прямого хода получили систему с матрицей верхнего треугольного вида.

     Преобразования прямого хода можно записать в следующем кратком виде:

~

~ ~ ~

~

Обратный ход. Находим неизвестные х₄, х₃, х₂, х₁ по формулам (12):

из 4-го уравнения находим: х₄ = –  » –11,1643,

из 3-го, 2-го и 1-го уравнений последовательно находим:

х₃ = 0,2572 – 1,5142 x₄ » 17,1622,

х₂ = –0,4286 + 1,8571 x₄ + x₃ » –3,9996,

х₁ = –0,5 + 2,5 x₄ + 0,5 x₃ – 1,5 x₂  » –13,8303.

Если округлить значения полученного решениядо 2-x десятичных знаков после запятой, получим: х₁»–13,83; х₂» –4,00; х₃ » 17,16;  х₄ » –11,16.

2). Для вычисления определителя матрицы системы по формуле (13) перемножаем ведущие элементы всех 3-х шагов на элемент a₄₄′′′ » –0,5144:

det A = a₁₁·a₂₂′·a₃₃′′·a₄₄′′′ = 2 (–3,5) (–5) (–0,5144) = – 18,004 » –18.

Ответы:

1) решение системы: х₁» –13,83; х₂» –4,00;  х₃ » 17,16;  х₄ » –11,16;

2) определитель матрицы системы det A » –18.

Задача 4. Дана таблица значений функции f (x):

Требуется:

   1) по табличным данным построить для функции f (x) интерполяционный полином 4-го порядка в форме Лагранжа и привести его к стандартному виду целого многочлена;

2) используя полученный полином, вычислить приближенное значение

функции f (x) в точке  = 0,2.

Решение.

1). Все значения x_i в данной таблице образуют равномерную сетку:

x_i = x₀ + ih, где x₀ = – 0,5, h = 0,3 для i = 0, 1, 2, 3, 4. Построим интерполяционный полином Лагранжа 4-го порядка по формуле (14), подставляя в нее значения x_i  и y_i = f (x_i):

      Приведем полученный полином к стандартному виду, произведя упрощения:

L₄(x) » 16,572(x⁴ – x³ + 0,15x² + 0,05x – 0,0056) – 110,016 (x⁴ – 0,7x³ – 0,21x² +

+ 0,167x – 0,014) + 152,691(x⁴ – 0,4x³ – 0,39x² + 0,086x + 0,028) – 31,537(x⁴ –

– 0,1x³ – 0,39x² – 0,031x + 0,007) + 12,8616 (x⁴ +0,2x³ – 0,21x² – 0,022x +

–  + 0,004) = 40,5716x⁴ + 5,0888x³  – 24,3618x² – 3,7180x + 5,5535

(промежуточные вычисления производим, используя 4 десятичных знака после запятой).

2). Вычислим значение f ( ) с округлением до 3-го знака после запятой:

f ( ) = f (0,2) » L₄(0,2) »

» 40,5716·0,2⁴ + 5,0888·0,2³  – 24,3618·0,2² – 3,7180·0,2 + 5,5535» 3,941.

Ответы:

1) интерполяционный полином в стандартном виде:

L₄(x) » 40,5716x⁴ + 5,0888x³  – 24,3618x² – 3,7180x + 5,5535;

2) значение f ( ) = f (0,2) » 3,941.

Задача 5. Дан определенный интеграл . Составить таблицу значений подынтегральной функции в точках  x_i = 1 + ih, где i  = 0, 1, …,10  с шагом

h = 0,1 и вычислить приближенное значение интеграла, используя эту таблицу и квадратурную формулу Симпсона.

Решение. Для построения таблицы значений подынтегральной функции  вычисляем значения x_i = 1 + ih, где i  = 0, 1, …,10, h = 0,1, затем значения функции в этих точках f (x_i) = e^3xi / (10 x_i) и заносим их в таблицу. Для удобства дальнейших вычислений значения функции во внутренних точках x_i с нечетными номерами i смещаем влево, а с четными номерами – вправо.

Вычислим определенный интеграл, используя формулу Симпсона. Для этого подставим в формулу (15)  2m = 10, h = 0,1, y_i = f (x_i) и получим:

» 0,0333·(2,0086 + 20,1714 + 4 (2,4648 + 3,8002 + 6,0011 + 9,6484 + 15,7299) + + 2·(3,0499 + 4,7633 + 7,5944 + 12,3004)) » 7,606.

Ответ: » 7,606.

Задача 6.Дана задача Коши для обыкновенного дифференциального уравнения: y′ = 3x – 2y², y(0) = 0,4. Решить задачу при помощи метода Рунге-Кутты на промежутке [0; 0,5] с шагом h = 0,1.

Решение. Для решения задачи потребуется сделать 5 шагов, т.к.

.

Результаты промежуточных расчетов по формулам (16) – (17)для i = 0, 1,..,4 будем вносить в рабочую таблицу. Описание расчетов 0-го шага приведем подробно.

      Номер шага i = 0, x₀ = 0, y(x₀) = y₀ = 0,4.

      Вычислим величины k₁, k₂, k₃, k₄:

;

;

;

;

следующее значение y₁ находим по формуле (17):

y₁  = y₀ + (k₁ + 2k₂ + 2k₃ + k₄) » 0,3846.

     Далее фиксируем номер шага i = 1, x₁= 0 +ih = 0,1, y(x₁) = y₁ » 0,3846, вычисляем величины k₁, k₂, k₃, k₄ по формулам:

; ;

;

 и значение y₂.

   Затем выполняем все это для i = 2, 3, 4, вычисляя на каждом шагу значения k₁, k₂, k₃, k₄ по формулам (17), а следующее значение y_i+1 по формуле (16).

     Результаты расчетов по шагам оформляем в виде рабочей таблицы:

     В итоге получаем таблицу приближенных значений решения задачи Коши y (x) на промежутке [0; 0,5] с шагом h = 0,1.

Ответ: приближенное решение задачи Коши получено в виде таблицы значений y_i » y(x_i).

Задача 7.Температура однородного стержня U = U (х, t) в сечении х в момент времени t  удовлетворяет уравнению теплопроводности. Используя метод сеток, найти решение смешанной задачи для уравнения теплопроводности  в области  при заданных условиях:

    начальное распределение температуры в стержне U (х, 0) = cosx;

    температура на концах стержня U (0; t) = 2 t + 1, U (1,5; t) = cos (1,5).

Решение. Для решения используем сетку с шагом h = 0,3 по переменной х и с шагом d = 0,015 по переменной t.Введем следующие обозначения.

Узлы сетки: x_i = ih для i = 0, 1, 2, 3, 4, 5, x_i [0; 1,5];

t_k = kd для k = 0,1, 2, 3, 4, 5, 6,  t_k [0; 0,09].

Сеточная функция: u_{i, k} = u(i, k) для i = 0, …5,  k = 0, …6, соответствующая приближенному решению задачи, т.е. u_{i, k} » U (x_i, t_k).

В данном случае шаг d = 0,015 удовлетворяет условию d = h²/6 = 0,09/6, поэтому расчетные формулы для вычисления всех значений функции U (x_i, t_k) = u_{i, k} в узлах сетки, покрывающей  область D, будем получать по формулам (21) – (24).

     Результаты вычислений будем заносить в таблицу, построенную в соответствии с заданной сеткой, и имеющей вид:

	i	0	1	2	……	5
k	xi tk	0	0,3	0,6	……	1,5
0	0	u0,0	u1,0	u2,0	……	u5,0
1	0,015	u0,1	u1,1	u2,1	……	u5,1
……	……	……	……	……	……	……
6	0,09	u0,6	u1,6	u2,6	……	u5,6

В этой таблице нужно заполнить затененные ячейки значениями u_{i, k}.

Порядок заполнения таблицы следующий.

1. Заполняем 0-й столбец (температура на левом конце стержня в момент времени t_k): u_{0, k} = 2 t_k + 1, для k = 0,1, 2, 3, 4, 5, 6.

2. Заполняем последний столбец (температура на правом конце стержня в момент времени t_k):  u_{5, k} = cos (1,5) для k = 0,1, 2, 3, 4, 5, 6.

3. Заполняем внутренние ячейки 0-й строки (начальная температура в точках с абсциссой x_i): вычисляем u_{i, 0} = cos (x_i) для i = 1, 2, 3, 4.

4. Вычисляем значения u_{i, k} во внутренних узлах сетки:

 для i = 1, 2, 3, 4, k = 0, 1, 2, 3, 4, 5.

Вычисления производятся по строкам, значения u_{i-1, k} , u_{i, k}  и u_{i+1, k}  берутся из предыдущей строки. Например, 1-я строка (k + 1 = 1):

u_{1, 1} = (u_{0, 0} + 4·u_{1, 0} + u_{2, 0})/6 = (1 + 4·0,9553 + 0,8253)/6 » 0,9411;

u_{2, 1} = (u_{1, 0} + 4·u_{2, 0} + u_{3, 0})/6 = (0,9553 + 4·0,8253 + 0,6216)/6 » 0,8131;

u_{3, 1} = (u_{2, 0} + 4·u_{3, 0} + u_{4, 0})/6 = (0,8253 + 4·0,6216 + 0,3624)/6 » 0,6124;

u_{3, 1} = (u_{3, 0} + 4·u_{4, 0} + u_{5, 0})/6 = (0,6216 + 4·0,3624 + 0,0707)/6 » 0,3570.

     После заполнения 1-й строки переходим к заполнению 2-й, затем 3-й и т.д., до 6-й строки.

В результате расчетов получаем таблицу:

	i	0	1	2	3	4	5
k	xi tk	0	0,3	0,6	0,9	1,2	1,5
0	0	1	0,9553	0,8253	0,6216	0,3624	0,0707
1	0,015	1,03	0,9411	0,8131	0,6124	0,3570	0,0707
2	0,03	1,06	0,9346	0,8009	0,6032	0,3518	0,0707
3	0,045	1,09	0,9332	0,7903	0,5943	0,3469	0,0707
4	0,06	1,12	0,9355	0,7814	0,5857	0,3421	0,0707
5	0,075	1,15	0,9406	0,7745	0,5777	0,3375	0,0707
6	0,09	1,18	0,9478	0,7694	0,5705	0,3331	0,0707

Ответ: приближенное решение задачи получено в виде таблицы значений

u_{i, k}  » U (x_i, t_k) в узлах сетки области D.

Рекомендуемая литература

1. Волков Е. А. Численные методы: учебное пособие для вузов / Е. А. Волков. – М.: Наука,1987 248 с.

2. Вержбицкий В. М. Основы численных методов: учебник для вузов / В. М. Вержбицкий. – М.: Высш. шк., 2002.– 840 с.

3. Данко, П. Е. Высшая математика в упражнениях и задачах: учебное пособие для втузов. В 2 ч. Ч.2 / П. Е. Данко, А. Г. Попов, Т. Я. Кожевникова.– М.: Оникс: Мир и образование, 2005.– 416 с.

4. Воробьева, Г. Н., Данилова, А. Н. Практикум по численным методам / Г. Н. Воробьева, А.Н. Данилова.– М.: Высш. школа, 1979.– 184 с.

5. Мостовская Л.Г. Вычислительная математика: учебно-методическое пособие ч.1 / Л. Г. Мостовская, А.-В.И.Середа, – Мурманск, 2001. – 86 с.