Исследование устойчивости и точности аппроксимации в численном методе

3.2 Аналитическое и графическое представление устойчивости в МКР. Число Куранта

Следующим шагом знакомства студентов предлагается показать устойчивость аналитически. Основным параметром, отвечающим за устойчивость решения, является число Куранта. Число Куранта - это безразмерная величина, характеризующая соотношение шагов по времени и пространству. Для аналитического представления устойчивости в численном методе решим задачу о распределении температуры в стержне, нагретом в начальный момент времени в середине, с течением времени.

Любое численное решение начинается с выбора уравнения, подходящего для решения той или иной задачи. Наш пример не является исключением. Для начала выберем уравнение, которое подойдёт для описания физического процесса переноса теплоты в теле. В нашей задаче таким уравнением будет дифференциальное уравнение теплопроводности Фурье без источника теплоты:

. (3.1)

Для численного решения аппроксимируем уравнение в виде явной разностной схемы по четырёхточечному шаблону явной схемы, соответствующей рисунку 3.1.:

. (3.2)

Рисунок 3.1 - Четырёхточечный шаблон явной разностной схемы

В ходе математических преобразований выражаем значение температуры на следующем временном шаге t(k+1, i) и замечаем, что оно зависит от значений температуры на предыдущих шагах по координате и, так называемого, числа Куранта - важного показателя, от которого полностью зависит устойчивость решения:

, (3.3)

где - число Куранта.

Выражение (3.3) показывает, что разностное уравнение на каждом шаге по времени связано c предыдущим слоем по координате. На основе полученного выражения и будем выполнять решение нашей задачи.

Решение будем проводить на компьютере при помощи вычислительной программы Mathcad. Данные для численного решения представлены на рисунке 3.2, а программное решение разностной схемы представлено на рисунке 3.3.

Рисунок 3.2 - Данные для численного решения

Рисунок 3.3 - Программное решение разностной схемы

Начальное распределение температуры вдоль расчётной области и решение для двух моментов времени показано в соответствии с рисунком 3.4 сплошной, пунктирной и штриховой линиями соответственно. Физически такое поведение вполне естественно - с течением времени тепло из более нагретой области переходит в менее нагретую, а зона изначально высокой температуры остывает. Стоит отметить, что данное решение получено при значении числа Куранта С = 0.4.

Рисунок 3.4 - Решение линейного уравнения теплопроводности при значении числа Куранта С=0.4

Чтобы явно показать устойчивость в МКР слегка изменим соотношение шагов по пространству и по времени, т.е. изменим число Куранта, произведя вычисления чуть-чуть большим шагом по времени. Оказывается, что незначительное увеличение шага до приводит совсем к другим результатам (значение числа Куранта в этом случае С = 1.6.) Результаты представлены на рисунке 3.5.

Рисунок 3.5 - Решение линейного уравнения теплопроводности при значении числа Куранта С=1.6

По сравнению с истинными температурными полями последний график, соответствующий рисунку 3.5, не отражает реальное охлаждение тела: появляется в алгоритме численного метода так называемая разболтка решения. Эта разболтка как раз и является результатом неустойчивости явной схемы сеточной функции для выбранного соотношения шагов по пространству и времени, т.е. числа Куранта. Таким образом показывается важность знакомства специалиста, использующего численный метод решения задач, с такой характеристикой, как устойчивость.

Данный уровень подготовки (аналитическое и графическое представление устойчивости) я предлагаю для студентов любого уровня подготовки (студенты-бакалавры и студенты-магистры), как вводный для знакомства с сутью явления устойчивости при численном решении.

3.3 Гармонический анализ устойчивости (Фурье)

Следующим этапом я предлагаю для студентов-магистров углубиться в природу возникновения разболтки решения, которую мы наблюдали при решении задачи с числом Куранта С =1,6. Выяснить, существует ли ограничение на выбор соотношения шагов по пространству и времени? То есть когда мы пытаемся улучшить нахождение величин температур в узлах шаблона, уменьшая шаги по пространству, как будет вести себя численный алгоритм явного метода. Для этого возьмем аппроксимированное явной схемой дифференциальное уравнение теплопроводности Фурье без источника теплоты:

. (3.4)

Перепишем эту схему в форме алгоритма для вычисления температур в последующие моменты времени:

. (3.5)

Для того, чтобы исследовать устойчивость явной схемы, необходимо выяснить, как происходить эволюция погрешности решения задачи при переходе от предыдущего слоя к последующему слою. Для этого предположим, что на каком-то этапе вычисления по времени у сеточной функции появилась погрешность относительно точного решения задачи :

. (3.6)

Подставим соотношение (3.6) в явное разностное уравнение (3.5) с учетом того факта, что значения температур при точном вычислении сократятся, тогда явная разностная схема дает следующее уравнение эволюции погрешности :

. (3.7)

Для того, чтобы найти относительное изменение погрешности в следующих этапах вычислений по времени, надо полученное выражение (3.7) разделить на , в результате получим:

(3.8)

Теперь определим динамику правой части равенства (3.8). Очевидно, что, если правая часть уравнения меньше 1, то явная схема устойчива, т.е. погрешности вычислений температур в узлах схемы со временем не возрастают. Если правая часть равенства больше 1, то схема неустойчива, так как погрешность при «движении» от слоя к слою возрастает. Для того, чтобы конкретно найти от каких параметров зависит устойчивость, зная из гармонического анализа Фурье, что любую функцию можно представить как гармоническую функцию, предположим, что величина погрешности имеет гармоническую зависимость от координат вычислительного шаблона с некоторым волновым числом :

, (3.9)

где - мнимая единица.

Подставляя выражение (3.9) в уравнение (3.8), получим:

(3.10)

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

· если , а скобка равна нулю, то равно единице;

· если , то при любом значении скобки меньше единицы, значит ошибки не растут, то есть явная схема является условно-устойчивой;

· если , а скобка равна 2, то по абсолютной величине больше единицы, то ошибки растут, т.е. это как раз условие неустойчивости явной разностной схемы, которое было показано графике.

Проверим выводы, к которым мы пришли, анализируя значение относительной погрешности в эволюционной задаче с послойным переходом. Для этого зададимся такой последовательностью рассуждений. Если t (x,t) - точное решение краевой задачи, а - точное решение конечно-разностных уравнений, то погрешность метода - разность . Важно выяснить, как изменяется при и фиксированных и , а также при , (когда ячейки сетки измельчаются) и фиксированном значении . Рассмотрим поведение явной разностной схемы при разных значениях сеточного числа Куранта. Допустим, что до (n-1) -го шага по времени расчёты велись точно, а на следующем (n-ом) временном слое в узле по какой-то причине возникла погрешность ; дальнейшие расчёты расчеты вновь будут точными. Исходя из этого, с зависимости от числа Куранта получим следующее распределение погрешности.

При С = 2 и учёте распределения погрешности в слоях значения погрешности на слое вычисляются по формулам:

(3.11)

(3.12)

(3.13)

. (3.14)

При С = и учёте распределения погрешности в слоях значения погрешности на слое вычисляются по формулам:

(3.15)

(3.16)

. (3.17)

Результаты расчётов для трёх временных слоёв представлены в таблице 3.1.

Таблица 3.1 - Результаты расчёта погрешностей явных схем

С	Шаг	Погрешность в слое
		j-3	j-2	j-1	j	j+1	j+2	j+3
1	2	3	4	5	6	7	8	9
2	n	0	0	0	?	0	0	0
	n+1	0	0	?	-?	?	0	0
	n+2	0	?	-2?	3?	-2?	?	0
	n+3	?	-3?	6?	-7?	6?	-3?	?
1/2	n	0	0	0	?	0	0	0
	n+1	0	0	0.5?	0	0.5?	0	0
	n+2	0	0.25?	0	0.5?	0	0.25?	0
	n+3	0.125?	0	0.375?	0	0.375?	0	0.125?
1/6	n	0	0	0	?	0	0	0
	n+1	0	0	0.176?	0.67?	0.17?	0	0
	n+2	0	0.01?	0.24?	0.5?	0.24?	0.01?	0
	n+3	0.01?	0.05?	0.245?	0.42?	0.245?	0.05?	0.001?

Для наглядности построим график распределения погрешности в слоях, соответствующий рисунку 3.6.

Рисунок 3.6 - График распределения погрешности в слоях в зависимости от значения числа Куранта

Видно, что распределение погрешности по узловым точкам сетки, а также значение местной ошибки существенно зависят от числа Куранта. Если при С=1/4 и С=1/6 начальная погрешность равномерно распределяется по всем узловым точкам сетки и тем быстрее, чем меньше значение числа Куранта, то при С = 2 ошибка Е быстро возрастает по абсолютной величине в соседних узловых точках. Изменение погрешности по знаку говорит о том, что изменение узловых значений температур во времени приобретает характер расходящихся колебаний, тогда как при и узловые температуры изменяются апериодически. Исходя из этого, можно утверждать, что выход изменения температур на режим колебаний с возрастающей амплитудой можно характеризовать, как потерю устойчивости численной схемы.

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

Говоря об устойчивости численной схемы, мы говорим об её нечувствительности к небольшим изменениям начальных данных, коэффициентам уравнений и случайным локальным погрешностям. Устойчивость нельзя связывать с ошибками округления, которые в любом случае будут неизбежны при выполнении вычислений; устойчивость - внутреннее свойство самой системы разностных уравнений.

Как мы уже отмечали, что конечная разностная схема устойчива, если для малых возмущений входных данных (начально-краевых условий и правых частей) конечная разносная схема обеспечивает малые возмущения сеточной функции , т.е. решение с помощью конечно-разностной схемы находится под контролем входных данных.

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

Из курса теплотехники и математической физики известно, что решение начально-краевых задач может быть представлено в виде суммы произведений двух функций:

, (3.18)

где л_n - это собственные значения.

- это собственные значения функции, получаемые из решения соответствующей задачи Штурма-Лиувиля. То есть решение может быть представлено в виде суперпозиции отдельных гармоник:

, (3.19)

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

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

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

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

, (3.20)

где - амплитудная часть.

Амплитудную часть представим в виде произведения двух функций, опираясь на гармонический анализ (Фурье):

, (3.21)

где - размерный и постоянный сомножитель амплитудной части;

k - показатель степени (соответствующий номеру временного слоя) сомножителя, зависящего от времени.

Тогда подставив (3.21) в конечно-разностную схему, можно по модулю оценить отношение амплитудных частей на соседних временных слоях.

Однако поскольку операция суммирования линейна и собственные функции ортогональны для различных индексов суммирования, то в конечно-разностную схему вместо сеточных значений достаточно подставить одну гармонику разложения (3.21) (при этом у амплитудной части убрать индекс n), т.е.:

, (3.22)

, (3.23)

, (3.24)

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

, (3.25)

т. е. условие (3.25) является необходимым условием устойчивости.

3.4 Спектральный анализ (Фурье)

Данный способ исследования устойчивости я предлагаю студентам-магистрам, которым необходимы углублённые, специальные знания по теме устойчивого решения. Он требует от студента глубоких знаний математического анализа, гармонического анализа (Фурье) и умения решать нетривиальные дифференциальные уравнения.

Рассмотрим разностную схему, записанную в каноническом виде. То есть у нас производная по времени представлена, как конечная разность:

, (3.26)

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

B - это оператор, от значения которого зависит явная или неявная схема перед нами: если оператор единичный, то мы получим явную разностную схему, если не единичный то неявную разностную схему.

Для явной схемы оператор L всегда может быть представлен, как:

, (3.27)

где - это неопределённые коэффициенты;

- принадлежит шаблону рассматриваемой разностной схемы.

Представим U, как функцию, от которой существует преобразование Фурье:

. (3.28)

Умножим выражение (3.28) на комплексную экспоненту и возьмем интеграл, получив тем самым преобразование Фурье от функции U:

. (3.29)

Возьмем преобразование Фурье от левой и правой части выражения (3.26) по пространственным координатам:

, (3.30)

где - это Фурье образ.

Стоит отметить, что в выражении (3.25) мы не берем Фурье-преобразование от f и, следовательно, не добавляем получившееся бы преобразование в выражение 3.30, которое превратило бы выражение (3.30) в неоднородное уравнение. Делаем мы так, потому что известно, что общее решение неоднородного уравнения есть сумма частного решения неоднородного уравнения и общего решения однородного решения. Все свойства решения в данном случае определяются свойствами решения однородного уравнения. В ходе преобразований получаем обычное разностной уравнение с постоянными коэффициентами. Данное разностное уравнение первого порядка теперь поддаётся решению. Известно, что для функции U с волной разностное решение будет иметь:

, (3.31)

где - это спектр оператора послойного перехода.

Спектр оператора послойного перехода находится, как корень уравнения:

, (3.32)

где k - это число Фурье, .

Разберёмся в чём практический смысл оператора послойного перехода лямбда. Для этого рассмотрим две близкие задачи. Постановка задач будет одинаковая в первом и втором случае в соответствии с выражением 3.33. Но предположим, что первую задачу мы будем решать на идеальном компьютере с бесконечной длиной мантиссы, т.е. без учёта ошибок округления, а вторая задача рассчитывается на реальном компьютере с ошибками округления.

. (3.33)

Введём погрешность решения:

. (3.34)

Также покажем, что в случае решения задачи на реальном компьютере у нас «набегает» и в значение некоторая погрешность:

. (3.35)

Отнимем от (3.33) выражение (3.35). Тем самым мы находим ошибки округления, которые возникли за счет перехода от идеальной вычислительной системы к реально работающему компьютеру с конечным числом знаков в мантиссе:

. (3.36)

Пусть у нас в начальных условиях присутствует небольшая случайная погрешность в виде функции от пространственной переменной в начальный момент времени. Применим преобразование Фурье к данной функции погрешности, тем самым переведя её из пространственной области в частотную область, за которую отвечает k из выражения (3.32). Мы получили спектр: каждому k у нас соответствует возмущение начальных данных. Каждому k соответствует , являющаяся функцией от k и i. Тогда - коэффициент усиления той или иной гармоники с волновым числом k.

Представим, что наша сеточная функция U в любой момент времени (или на каждом слое по времени) представляет собой периодическую функцию. Тогда после того, как мы нашли преобразование Фурье, мы понимаем, что спектр U частотной области - это конечная сумма гармоник. И количество гармоник, участвующих в начальном возмущении - конечное число.

Исходя из этого, можно предположить, что если каждая из гармоник будет затухать, то и их сумма будет затухать. Возмущения не будут накапливаться. Потребуем, чтобы выполнялось условие:

. (3.37)

Тогда ни одна гармоника возрастать не будет. Следовательно, не будет возрастать и их сумма.

Если задаться вопросом, как влияют ошибки округления, то следует добавить к выражению (3.37) дополнительное слагаемое . Величина - некая величина, зависящая от разрядности компьютера, на котором происходят вычисления и от того, какая у нас правая часть в выражении (3.36). В следствие мы получим новое уравнение (3.38) и подойдём к формулировке спектрального признака устойчивости Фон-Неймана или, как его еще называют, признака Фурье:

. (3.38)

Признак устойчивости Фурье: если мы находим спектр оператора послойного перехода (коэффициент усиления перед каждой гармоникой Фурье) и если у нас для любого k выполнено неравенство (3.38), то наша схема устойчива. Она удовлетворяет условию устойчивости.

Также из признака устойчивости Фурье следует определение устойчивости. Заметим, что - величина маленькая и не зависит от сеточных параметров, h, и , а зависит только от свойств правой части выражения (3.36) и разрядности вычислительной машины. Возведём в выражении (3.37) в n-ую степень, где n - это текущее время, делённое на шаг по времени:

. (3.39)

Перейдём в пределе при , чтобы показать, что схема выдерживает последовательность вычислений на всё более и более мелких сетках. Получим выражение:

. (3.40)

Из выражения (3.40) следует, что любая погрешность (связанная с округлением в том числе) не будет возрастать до бесконечности (она будет всегда ограничена). Из этого положения и будет следовать устойчивость. Таким образом мы обосновали спектральный признак Фурье.

Следующим пунктом заметим, что, при фиксированном шаге h в выражении (3.32), произведение - это параметр. Значит мы ищем решение однородного уравнения, которое описывает эволюцию погрешности, связанную с округлением и наличием ошибок в начальных данных. Представим решение такого уравнения в виде уравнения:

, (3.41)

где m - номер точки;

.

Так как, то тоже может быть любым, но так как на сетке возможны не все гармоники, поэтому .

Если мы решение в виде (3.41) подставим в нашу разностную схему, то мы сумеем получить простое алгебраическое уравнение для спектра оператора послойного перехода. Применим спектральный признак устойчивости к исследованию разностных схем.

Например, явная разностная схема для уравнения теплопроводности имеет вид:

. (3.42)

Подставляем решение в виде уравнения (3.31), выносим произведение за скобку и получаем:

, (3.43)

где .

Отсюда очевидно находится спектр оператора послойного перехода. Разделим на ненужное и умножим левую и правую части на получаем:

, (3.44)

где .

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

. (3.45)

Следующим вариантом применения признака устойчивости рассмотрим анализ неявной схемы, когда производную по пространственной переменной мы берем на «верхнем» слое:

. (3.46)

Уравнение на спектр будет получаться таким же, только будет в другом месте:

. (3.47)

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

. (3.48)

Выражение (3.48) обладает совсем другим свойствами, т.к. число Фурье (k) естественно положительная величина. Соответственно при любом значении значение по модулю не превосходит единицы.

В случае, когда , у нас получается величина строго меньше единицы. Следовательно, все гармоники затухают. В этом случае по спектральному признаку наша разностная схема для уравнения теплопроводности получается безусловно устойчива, т.е. у нас нет никаких строгих условий для соотношения шагов по времени и пространству. Разностная схема условно устойчива.

3.5 Схема Кранка Николсона

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

, (3.49)

где - число Куранта.

Если записать уравнение (3.49) через сеточное число Фурье , оно примет вид:

. (3.50)

Как известно, данное уравнение обладает тем достоинством, что решение на верхнем временном слое t_k+l получается сразу (без решения СЛАУ) по значениям сеточной функции на нижнем временном слое t^k, где решение известно (при k = 0 значения сеточной функции формируются из начального условия). Но эта же схема обладает существенным недостатком, поскольку она является условно устойчивой. С другой стороны, неявная конечно-разностная схема, записанная форме выражения

, (3.51)

приводит к необходимости решать систему линейных алгебраических уравнений (СЛАУ), но зато эта схема абсолютно устойчива.

Проанализируем эти схемы (3.50) и (3.51). Пусть точное решение задачи возрастает по времени при , которое пока неизвестно. Тогда, в соответствии с явной схемой (3.50), разностное решение будет заниженным по сравнению с точным. Это происходит из-за того, что температура в следующий момент времени определяется по меньшим значениям сеточной функции на предыдущем временном слое, поскольку решение является возрастающим по времени. Для неявной схемы (3.51) на возрастающем решении, наоборот, решение завышено по сравнению с точным. Значение температуры определяется по значениям сеточной функции на верхнем временном слое.

Если рассмотреть изменения температур на убывающем участке, то картина изменяется противоположным образом: явная конечно-разностная схема завышает решения, а неявная занижает, в соответствии с рисунком 3.7.

Рисунок 3.7 - Графики, условно показывающие решение задачи явной схемой, решение неявной схемой и точное решение

Исторически на основе этого анализа возникла идея совместить преимущества этих двух схем, построив более точную совместную неявно-явную конечноразностную схему с весами при пространственных конечно-разностных операторах. Причем при измельчении шагов точное (неизвестное) решение может быть взято в «вилку» сколь угодно узкую, так как если явная и неявная схемы аппроксимируют дифференциальную задачу и эти схемы устойчивы, то при стремлении сеточных характеристик ф и к нулю решения по явной и неявной схемам стремятся к точному решению с разных сторон. Проведенный анализ дал блестящий пример так называемых двусторонних методов, исследованных В. К. Саульевым.

Рассмотрим неявно-явную схему с весами для простейшего уравнения теплопроводности:

, (3.52)

где и - вес неявной части конечно-разностной схемы;

и-1 - вес для явной части.

Причем . При и=1 имеем полностью неявную схему, при и = 0 - полностью явную схему, а при и=1/2 - схему Кранка-Николсона.

В соответствии с гармоническим анализом, для схемы (3.52) получаем неравенство:

, (3.53)

из которого получаем:

, (3.54)

где правое неравенство выполнено всегда.

Левое неравенство имеет место для любых значений , если . Если же вес и лежит в пределах , то между и и из левого неравенства устанавливается связь:

(3.55)

При , (3.56) являющаяся условием устойчивости неявно-явной схемы с весами (3.52), когда вес находится в пределах .

Таким образом, неявно-явная схема с весами абсолютно устойчива при и условно устойчива с условием (3.44) при .

Рассмотрим порядок аппроксимации неявно-явной схемы с весами, для чего разложим в ряд Тейлора в окрестности узла () на точном решении значения сеточных функций по переменной , а , по переменной х и полученные разложения подставим в (3.52). Получим выражение:

. (3.57)

В этом выражении дифференциальный оператор от квадратной скобки в соответствии с дифференциальным уравнением равен дифференциальному оператору , в соответствии с чем вышеприведенное равенство приобретает вид:

. (3.58)

После упрощения получаем:

, (3.59)

откуда видно, что для схемы Кранка-Николсона (и = 1/2) порядок аппроксимации схемы (3.42) составляет , т.е. на один порядок по времени выше, чем для обычных явных или неявных схем. Таким образом, схема Кранка-Николсона при абсолютно устойчива и имеет второй порядок аппроксимации по времени и пространственной переменной х.

Используем в уравнение (3.52) подстановку . Но в то же время его нужно решить для трех "еще не вычисленных" значений и . Это возможно, если все значения перенести в левую часть уравнения. Затем упорядочим члены уравнения (3.52) и в результате получим неявную разностную формулу:

. (3.50)

Члены в правой части формулы (3.50) известны. Таким образом, формула (3.50) имеет вид линейной трехдиагональной системы . Шесть точек, используемых в формуле Кранка-Николсона (3.50), вместе с промежуточной точкой решетки, на которой основаны численные приближения, в соответствии с рисунком 3.8.

Рисунок 3.8 - Шаблон (разностная схема) метода Кранка-Николсона

Иногда в формуле (3.60) используется значение . В этом случае приращение по оси t равно , тогда формула (3.60) упрощается и принимает вид:

. (3.61)

Граничные условия используются в первом и последнем уравнениях, то есть в

(3.62)

(3.63)

Уравнения (3.62) и (3.63) особенно привлекательны при записи в форме трехдиагональной матрицы .

Если метод Кранка-Николсона реализуется на компьютере, то линейную систему можно решить либо прямым методом, либо итерационным.

Пример. Найдём с помощью метода Кранка-Николсона решение уравнения:

, (3.64)

где ;

.

Уравнение (3.64) имеет начальное условие:

(3.65)

где ;

.

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

, (3.66)

. (3.67)

Решение будем искать в системе автоматизированного проектирования Mathcad. Для простоты возьмем шаг Дх = h = 0,1 и Дt = к = 0,01. Таким образом, соотношение r =1. Пусть решетка имеет n=11 столбцов в ширину и m=11 рядов в высоту.

Решения для данных параметров отразим в таблице 3.2.

Таблица 3.2 - Значения , полученные методом Кранка-Николсона

xi	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
ti
0	0	1.1180	1.5388	1.1180	0.3633	0	0.3633	1.1180	1.5388	1.1180	0
0.01	0	0.6169	0.9288	0.8621	0.6177	0.4905	0.6177	0.8621	0.9288	0.6169	0
0.02	0	0.3942	0.6480	0.7186	0.6800	0.6488	0.6800	0.7186	0.6480	0.3942	0
0.03	0	0.2887	0.5067	0.6253	0.6665	0.6733	0.6665	0.6253	0.5067	0.2887	0
0.04	0	0.2331	0.4258	0.5560	0.6251	0.6458	0.6251	0.5560	0.4258	0.2331	0
0.05	0	0.1995	0.3720	0.4996	0.5754	0.6002	0.5754	0.4996	0.3720	0.1995	0
0.06	0	0.1759	0.3315	0.4511	0.5253	0.5504	0.5253	0.4511	0.3315	0.1759	0
0.07	0	0.1574	0.2981	0.4082	0.4778	0.5015	0.4778	0.4082	0.2981	0.1574	0
0.08	0	0.1419	0.2693	0.3698	0.4338	0.4558	0.4338	0.3698	0.2697	0.1419	0
0.09	0	0.183	0.2437	0.3351	0.3936	0.4137	0.3936	0.3351	0.2437	0.1283	0
0.1	0	0.1161	0.2208	0.3038	0.3570	0.3753	0.3570	0.3038	0.2208	0.1161	0

Величины, полученные методом Кранка-Николсона, достаточно близки к аналитическому решению:

, (3.68)

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

Таблица 3.2 - Точные значения при t = 0.1

xi	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
t11
0.1	0	0.1153	0.2192	0.3016	0.3544	0.3726	0.3544	0.3016	0.2192	0.1153	0

Рисунок 3.9 - Численное решение задачи (3.54) в трёхмерном виде

Проанализируем полученные численным методом значения температуры и сравним их с истинными значениями параметра. Для этого построим график, на котором представим значения аналитического решения и численного, соответствующий рисунку 3.10.

Рисунок 3.10 - Аналитическое и численное решение задачи (3.54)

Исходя из сравнения двух графиков, видно, что, используя схему Кранка-Николсона, мы получаем значения параметра, отличающиеся от истинных в среднем на 0,7%. Данная погрешность решения более чем удовлетворяет требованиям точности при широком спектре задач современной строительной отрасли.

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

ЗАКЛЮЧЕНИЕ

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

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

· найден материал для знакомства студентов с характеристиками численных методов, показана их значимость при решении задач численным методом;

· показан способ аппроксимации дифференциальных уравнений конечно-разностными - схемами;

· показан способ определения порядка аппроксимации: разложение функции в ряд Тейлора;

· показаны способы увеличения порядка аппроксимации: метод фиктивных узлов, метод уменьшения невязки, метод неопределённых коэффициентов;

· показаны формы, в которых проявляется устойчивость при численном решении: графическое представление устойчивого и неустойчивого численных решений, аналитическое представление на примере числа Куранта, гармонический анализ устойчивости (Фурье), спектральный анализ (Фурье);

· описана природа устройства существования устойчивости в конечно-разностных схемах;

· показаны способы воздействия на устойчивость численного решения для получения корректного решения задачи;

· произведено распределение информации по характеристикам численных методов в зависимости от уровня подготовки студентов;

· предложены этапы обучения студентов анализу устойчивости конечно-разностных схем.

Дальнейшие перспективы развития исследования данной темы я вижу в разработке заданий для студентов разных уровней подготовки (студенты-бакалавры и студенты-магистры) с целью приобрести умение решать дифференциальные уравнения в численном виде и использовать характеристики численных методов для получения корректного решения.

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1. Телегин А. С. Тепломассоперенос / А. С. Телегин, В. С. Швыдкий, Ю. Г. Ярошенко. - Москва: Металлургия, 1995. - 400 с.

2. Лекция 6: Устойчивость разностных схем для эволюционных уравнений и решение [Электронный ресурс]: видеохостинг youtube.com. - Режим доступа: https://www.youtube.com/watch?v=2CL5FbFKFtk&t=8s.

3. Пискунов Н. С. Дифференциальное и интегральное исчисления: учеб. пособие / Н. С. Пискунов. - Москва: Наука, 1978. - 456 с.

4. Кирьянов Д. В. Mathcad 12: учеб. пособие / Д. В. Кирьянов. - Санкт-Петербург: БХВ-Петербург, 2005. - 576 с.

5. Колдаев В. Д. Численные методы и программирование: учеб. пособие / Д. В. Колдаев. - Москва: Форум, 2009. - 336 с.

6. Нагрузки на ролики при вибрационном смещении опор на МНЛЗ/ А.И. Буркова, Н.Г. Баширов, В.А. Куницкий // Инновационные технологии в науке и образовании - 2018. - №1. - С.133-136

7. Обучение студентов анализу устойчивости численных методов/ Н.Г. Баширов, А.И. Буркова, И.А. Баширова // Вузовская наука региону: материалы шестнадцатой всероссийской научно-технической конференции. - Вологда: ВоГУ, 2018.

8. Куницкий В.А. Обучение студентов анализу устойчивости численных методов В.А. Куницкий, Молодые исследователи - регионам: материалы международной научной конференции. / Мин-во обр. и науки РФ; Вологод. гос. ун-т. - Вологда: ВоГУ, 2018 (принято в печать).

Размещено на Allbest.ru

...