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

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

1. ТЕОРЕТИЧЕСКИЕ ОСНОВЫ ЧИСЛЕННЫХ МЕТОДОВ. ЦЕЛИ И ЗАДАЧИ ДИПЛОМНОЙ РАБОТЫ

1.1 Виды численных методов

1.2 Характеристики численных методов

1.3 Точность аппроксимации (порядок аппроксимации)

1.4 Устойчивость

1.5 Цели и задачи

2. ИССЛЕДОВАНИЕ ТОЧНОСТИ АППРОКСИМАЦИИ В ЧИСЛЕННЫХ МЕТОДАХ

2.1 Аппроксимация дифференциальных уравнений. Способ определения порядка аппроксимации в разностных схемах

2.2 Способ фиктивных узлов

2.3 Метод уменьшения невязки

2.4 Метод неопределённых коэффициентов, как способ повышения порядка аппроксимации

3. ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ В ЧИСЛЕННЫХ МЕТОДАХ

3.1 Знакомство с устойчивостью в МКР

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

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

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

3.5 Неявно-явная схема Кранка-Николсона

ЗАКЛЮЧЕНИЕ

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

ВВЕДЕНИЕ

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

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

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

В подобных случаях пользуются численными методами. Численные методы - математические методы приближённого решения математических задач, основанные на переходе от непрерывных функций к дискретным и аппроксимации (замене) исходных дифференциальных уравнений алгебраическими. Подобные преобразования позволяют представить область изменения параметров общностью (разностной сеткой) взаимосвязанных точек (узлов сетки) со значениями параметра в каждой из них и записать формулы для определения значения параметра в любой точке. Значение параметра (функции) при такой замене будет зависеть от номера слоя и от значений параметров на соседних слоях. Общность слоёв образуют разностную сетку - дискретный аналог представления времени и пространства. Разностная сетка включает в себя множество узлов у каждого из которых есть свой номер слоя по времени и координате. А расстояние между узлами - шаг сетки, значение которого для времени и пространства может быть разным.

Я прохожу обучение на кафедре «Теплогазоснабжения и вентиляции». В учебном плане для групп этой кафедры СТ-31 и СТЭ-31 есть курс «Численные методы». Для успешного его освоения от студентов требуется знание основ математики и высокий уровня компьютерной грамотности.

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

В связи с этим возникла проблема эффективного и экономичного по времени обучения студентов анализу данных характеристик.

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

1. ТЕОРЕТИЧЕСКИЕ ОСНОВЫ ЧИСЛЕННЫХ МЕТОДОВ. ЦЕЛИ И ЗАДАЧИ ДИПЛОМНОЙ РАБОТЫ

1.1 Виды численных методов

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

· точность аппроксимации;

· устойчивость;

· сходимость.

Принято выделять два основных вида численных методов: метод конечных разностей (МКР) и метод конечных элементов (МКЭ). У каждого из них есть свои преимущества и недостатки.

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

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

· на плоскости в рассматриваемой области А в которой ищется решение, строится сеточная область А (i, k) в соответствии с рисунком 1.1, состоящая из ячеек размером (i, k) (I, k - шаг сетки) и являющаяся приближённым изображением данной области А;

· осуществляется переход от непрерывной функции-решения рассматриваемого дифференциального уравнения к дискретной функции, т.е. дифференциальное уравнение аппроксимируется (заменяется) в узлах сетки А (i, k) конечно-разностным уравнением;

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

· краевая задача заменяется дискретной краевой задачей (разностной схемой), т.е. системой алгебраических уравнений;

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

Рисунок 1.1 - Аппроксимация объекта сеточной областью

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

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

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

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

· МКР позволяет получить расчетные значения во всех узлах области, а не только интегральные значения по поперечному сечению;

· построение разностной схемы происходит быстро для простых задач.

Недостатками МКР являются:

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

· расчеты невозможно выполнить вручную; для их осуществления необходим ЭВМ с рабочим алгоритмом вычисления;

· сильная «чувствительность» к аппроксимации исходных данных и затруднение при построении разностной схем для сложных геометрических тел.

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

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

· задачи гидромеханики;

· задачи механики и прочности.

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

Рассмотрим построение дискретной модели на примере одномерной задачи о распределении температуры в стержне в соответствии с рисунком 1.2.

Рисунок 1.2 - График функции, описывающей распределения температуры в стержне

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

· в рассматриваемой области фиксируется конечное число точек, эти точки называются узловыми или просто узлами в соответствии с рисунком 1.3;

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Рисунок 1.3 - Узловые точки

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

· область определения непрерывной величины разбивается на конечное число подобластей, называемых элементами (или конечными элементами), изображёнными на рисунке 1.4; эти элементы имеют общие узловые точки и в совокупности аппроксимируют форму области; разбиение области на элементы может быть произведено двумя различными способами:

а) ограничение каждого элемента двумя соседними узловыми точками с образованием четырёх элементов в соответствии с рисунком 1.4а;

б) разбиение области на два элемента, каждый из которых содержит три узла в соответствии с рисунком 1.4б;

Рисунок 1.4 - Разбиение непрерывной функции на конечное число подобластей

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

При построении дискретной модели непрерывной величины, определённой в двух или трёхмерной области, основная концепция МКЭ используется аналогично. В двумерном случае элементы описываются функциями, зависящими от переменным по осям абсцисс и ординат. При этом чаще всего рассматриваются элементы в форме треугольника или четырёхугольника. Функции элементов изображаются плоскими поверхностями в соответствии с рисунком 1.5 или криволинейными поверхностями в соответствии с рисунком 1.6.

Рисунок 1.5 - Изображение функции плоскими поверхностями.

Рисунок 1.6 - Изображение функции криволинейными поверхностями.

устойчивость аппроксимация численный

Функция элемента будет представляться плоскостью, если для данного элемента взято минимальное количество узловых точек, которое для треугольного элемента равняется трём, а для четырёхугольного - четырём. Если используемое число узлов больше минимального, то функции элемента будет соответствовать криволинейная поверхность. Кроме того, избыточное число узлов позволяет рассматривать элементы с криволинейными границами. Окончательной аппроксимацией двумерной непрерывной величины T (x, y) будет служить совокупность кусочно-непрерывных поверхностей, каждая из которых определяется на отдельном элементе с помощью значений T (x, y) в соответствующих узловых точках. Искомые узловые значения Т(х) должны быть «отрегулированы» таким образом, чтобы обеспечивалось «наилучшее» приближение к истинному распределению температуры. Это «регулирование» осуществляется путём минимизации некоторой величины, связанной с физической сущностью задачи. Если рассматривается задача распространения тепла, то минимизируется функционал, связанный с уравнением теплопроводности.

Преимущества метода конечных элементов:

· свойства материалов смежных элементов не должны быть обязательно одинаковыми; это позволяет применять метод к телам, составленным из нескольких материалов;

· криволинейная область может быть аппроксимирована с помощью прямолинейных элементов или описана точно с помощью криволинейных элементов; таким образом, методом можно пользоваться не только для областей с «хорошей» формой границы;

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

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

Тем не менее МКЭ тоже не свободен от недостатков:

· время, необходимое для расчетов, а также требования к аппаратным средствам компьютера и объему носителей информации в несколько раз превышают аналогичные требования для МКР; для решения задач этим методом требуется как минимум высокопроизводительный 16-разрядный или 32-разрядный ПК; за редким исключением, применение программ, реализующих МКЭ, ограничивается плоскими задачами;

· поскольку геометрия канала, а также начальные и граничные условия задаются пользователем самостоятельно, время, необходимое для расчета, существенно больше, чем для МКР, где эти параметры более или менее фиксированы;

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

1.2 Характеристики численных методов

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

Эффективность применения ЭВМ во многом зависит от опыта, профессиональной квалификации и компьютерной грамотности специалиста. Одной из составляющих компьютерной грамотности на современном этапе является умение формализовать свои профессиональные знания и довести их до практического алгоритма использования. Это умение построить математическую модель технического процесса или изучаемого объекта, знание, какими численными методами может быть решена та или иная задача и умение выбрать наиболее рациональный из них и оценить достоверность полученных результатов. Также на некоторых этапах вычислительного эксперимента возникают погрешности, искажающие результаты вычислений. Поэтому оценка степени достоверности получаемых результатов в процессе вычислительных работ является важным вопросом. Рассмотрим последний пункт подробнее.

Исследуя источники ошибок (погрешностей) при использовании численных расчетов отметим ряд возможных причин некорректного решения:

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

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

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

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

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

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

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

Корректность - это, когда решение задачи существует, оно единственно и устойчиво относительно исходных данных из некоторого класса ее решений.

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

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

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

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

1.3 Точность аппроксимации (порядок аппроксимации)

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

Рассмотрим классический алгоритм аппроксимации краевой задачи. Область непрерывного изменения аргументов искомой функции t заменяется некоторым конечным множеством точек в этой области. Это множество называется разностной сеткой, сами точки - узлами сетки, а функции, определенные тем или иным способом на этой сетке, - сеточными функциями. Расположение узлов сетки в области может быть произвольным, оно определяется спецификой решаемой задачи. Поэтому можно выделить равномерные и неравномерные сетки. В сетках первого типа все узлы по одной координате расположены на равном расстоянии один от другого, называемым шагом сетки. Пример такой сетки для дифференциального уравнения теплопроводности показан на рисунке 1.7.

Рисунок 1.7 - Разностная стека для дифференциального уравнения теплопроводности

Здесь координаты любой точки на плоскости рисунка (х, т) определяются соотношениями:

, (1.1)

, (1.2)

где j = 1, 2, …;

- величина шага по пространству;

j - число шагов по толщине.

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

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

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

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

· метод неопределенных коэффициентов;

· метод интегральных тождеств (интегрально - интерполяционный).

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

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

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

(1.3)

на равномерной сетке:

x_i = i · h, i = 0, 1,..., N; t_j = j · ф j = 0, 1,..., M. M· ф = 1. (1.4)

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

(1.5)

Графическое представление сетки (1.5) показано на рисунке 1.8.

Рисунок 1.8 - Сетка аппроксимации с использованием фиктивных узлов

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

Рисунок 1.9 - Разностная схема

На сетке (1.4) разностная схема, представленная в соответствии с рисунком 1.9, обладает первым порядком погрешности аппроксимации по h и вторым по ф, а на сетке (1.5) -- вторым порядком погрешности аппроксимации по h и ф.

Решим задачи численно и проанализируем величину погрешностей, полученных на сетке (1.4) и сетке (1.5) между собой. Значения погрешности на сетке (1.4) и (1.5) представлены в соответствии с рисунком 1.10 и рисунком 1.11 соответственно.

Рисунок 1.10 - Погрешность решения задачи на сетке (1.4)

Рисунок 1.11 - Погрешность решения задачи на сетке (1.5)

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

1.4 Устойчивость

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

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

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

. (1.6)

Аппроксимацию будем вести четырехточечным шаблоном явной схемы без источника теплоты в соответствии с рисунком 1.12.

Рисунок 1.12 - Четырехточечная явная разностная схема

Получаем уравнение явной разностной схемы:

. (1.7)

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

(1.8)

или

, (1.9)

где параметр - называется числом Куранта.

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

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

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

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

Результаты численного решения с числом Куранта равным С = 1.6 представлены на рисунке 1.14.

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

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

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

1.5 Цели и задачи дипломной работы

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

· сбор и обработка уже имеющихся данных по характеристикам численных методов;

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

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

Для достижения этих целей подобраны задачи дипломной работы:

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

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

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

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

· познакомить с понятием аппроксимации;

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

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

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

2. ИССЛЕДОВАНИЕ ТОЧНОСТИ АППРОКСИМАЦИИ В ЧИСЛЕННЫХ МЕТОДАХ

2.1 Аппроксимация дифференциальных уравнений. Способ определения порядка аппроксимации в разностных схемах

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

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

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

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

Для примера возьмем дифференциальное уравнение теплопроводности:

. (2.1)

Область, в которой будет искаться решение аппроксимирована в виде разностной сетки, представленной на рисунке 2.1.

Рисунок 2.1 - Разностная стека для дифференциального уравнения теплопроводности

Здесь координаты любой точки на плоскости рисунка (х, т) определяются соотношениями:

, (2.2)

, (2.3)

где j = 1, 2, …;

- величина шага по пространству;

j - количество шагов по толщине.

Запишем следующее разложение для уравнения (2.1) по формуле Тейлора:

(2.4)

Запишем значение сеточной функции t в точке через нижний индекс j. Тогда перепишем выражение (2.4) в виде:

. (2.5)

Выразим первую производную:

. (2.6)

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

. (2.7)

Из сопоставления (2.7) с (2.6) определим интересующую нас погрешность произведённой аппроксимации:

, (2.8)

где .

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

Из выражения (2.8) невозможно получить точную величину погрешности, так как фактическое значение неизвестно, однако:

. (2.9)

Аналогичным образом, используя разложение по формуле Тейлора, получаем:

. (2.10)

Выражаем первую производную из уравнения (2.10):

, (2.11)

где .

Таким образом мы приходим к аппроксимации разностью назад (левой разностью):

. (2.12)

Погрешность данной аппроксимации снова имеет порядок и подчиняется условию:

. (2.13)

Важно подчеркнуть следующее: аппроксимации (2.7) и (2.12) основаны на разложении функции t в ряде Тейлора с точностью до линейных слагаемых. Благодаря этому они будут строго справедливы для линейных функций.

Вычитая (2.10) из (2.4), получаем:

. (2.14)

Данное действие позволяет вести нам аппроксимацию центральной разностью:

. (2.15)

Погрешность аппроксимации в данном случае будет:

. (2.16)

Так как порядок погрешности стал , то данное представление производной должно быть более точным, чем аппроксимации назад и вперёд. Данное утверждение становится очевидным, если отметить, что аппроксимация центральной разностью предполагает квадратичное изменение функции t между узлами и . К слову, стоит отметить, что в общем случае может возникнуть необходимость обеспечения «гладкой» стыковки отрезков шаблона параболы в процессе численного решения. Специалисту необходимо помнить этот нюанс при работе с подобной аппроксимацией.

Складывая разложения в ряд Тейлора (2.4) и (2.10), находим, что слагаемые с первой и третьей производными взаимно уничтожаются. В результате получаем формулу для нахождения второй производной:

. (2.17)

Исходя из этого, вторую производную можно аппроксимировать выражением:

. (2.18)

Погрешность аппроксимации (2.18) имеет порядок аппроксимации причем значение погрешности соответствует неравенству:

. (2.19)

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

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

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

.(2.20)

Рисунок 2.2 - Типичные шаблоны конечно-разностной схемы

Используя индексные обозначения и пренебрегая e(t) перезапишем выражение (2.20):

. (2.21)

где - сеточное число Фурье, составляющее половину от числа Куранта.

Для момента (n=1) распределение температур по сечению нагреваемого тела известно из начальных условий. Из-за этого расчёт температур на следующем временно слое не вызывает сложностей. Начиная «прогонять» расчёт для слоёв n=3,4,… и т.д., находим изменение всего температурного поля во времени. Так как новые температуры определяются из (2.20) явным образом, то данная разностная схема называется явной. То есть расчёт неизвестного значения параметра в узле на слое ведётся по известным значениям параметра в определённых узлах на предыдущем слое.

Используем шаблон, соответствующий рисунку 2.2б и запишем:

(2.22)

Используя индексные обозначения и пренебрегая e(t) перезапишем выражение (2.20):

. (2.23)

Легко заметить, что в этом случае температуры нового временного слоя явным образом не определяются. Здесь необходимо записать уравнение (2.21) для каждого узла сетки и затем решать совместно получившуюся систему уравнений. По этой причине разностная схема (2.23) называется неявной.

Совместим уравнения (2.21) и (2.23). В результате получим однопараметрическое уравнение конечно-разностных уравнений, определяемых на шеститочечном шаблоне, соответствующее рисунку 2.1.2в:

. (2.24)

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

и - числители разностных аппроксимаций второй производной .

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

Заменим производную по времени центральной производной и запишем:

. (2.25)

Есть множество различных схем для решения уравнения теплопроводности. Большинство из них можно получить из (2.24), просто варьируя значение .

Для оценки порядка разностной аппроксимации воспользуемся соотношением (2.4) и уравнениями:

(2.26)

и

. (2.27)

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

. (2.28)

Видно, что при формула (2.24) аппроксимирует уравнение теплопроводности с тем же порядком, что и явная схема. При порядок аппроксимации по времени увеличивается до двух, а для порядок аппроксимации по координате увеличивается до четырёх. В частном случае, когда , погрешность аппроксимации уменьшается до .

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

2.2 Способ фиктивных узлов, как метод повышения порядка аппроксимации

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

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

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

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

Рассмотрим вначале линейную задачу с постоянными коэффициентами:

(2.29)

и соответствующую ей разностную задачу:

. (2.30)

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

. (2.31)

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

Часто в данной ситуации пользуются способом создания фиктивного узла (ячейки), с помощью которого расширяют сеточную область. Подробнее рассмотрим данный способ. Введем в рассмотрение узлы с индексами 0 и N + 1, находящиеся формально за пределами рассматриваемой области. Пусть в них тоже определено значение сеточной функции. Тогда узел с индексом 0 выступает как внутренний узел, и в нем может быть записано соотношение:

. (2.32)

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

. (2.33)

Из этих двух уравнений осталось только исключить значение в фиктивном узле (фиктивной ячейке).

Также рассмотрим применение данного способа для уравнения разностного уравнения теплопроводности с невязкой, которая имеет первый порядок аппроксимации:

. (2.34)

Рассмотрим алгоритм повышения порядка аппроксимации граничных условий с первого до второго порядка .

Аналогично вводим вне рассматриваемого тела фиктивный узел и предполагаем, что уравнение теплопроводности справедливо для точек . Тогда разностное уравнение для явной схемы имеет вид:

. (2.35)

При значении можно записать:

. (2.36)

Следующим шагом для левого граничного условия аппроксимируем производную с помощью центральной разности, включая фиктивный узел:

. (2.37)

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

. (2.38)

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

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

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

2.3 Метод уменьшения невязки

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

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

Также, как в предыдущем пункте, разложим функцию температуры при помощи формулы Тейлора:

. (2.39)

Допустим, что у нас есть аппроксимированное граничное условие:

. (2.40)

А из уравнения теплопроводности получим выражение:

. (2.41)

Подставим уравнения (2.40) и (2.41) в формулу Тейлора и найдём:

+…. (2.42)

Заменяя производную по времени на границе и переходя к индексным обозначениям, снова получаем формулу:

. (2.43)

Поставим перед собой цель повысить порядок аппроксимации граничных условий третьего рода. Задаем граничные условия III рода на правой поверхности пластины:

, (2.44)

где - коэффициент теплоотдачи, ;

- температура среды, .

Используя разложение по формуле Тейлора получим:

. (2.45)

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

, (2.46)

где - сеточное число Био;

- сеточное число Фурье.

Порядок аппроксимации в данном случае составляет , что и является достигнутой нами целью по повышению порядка аппроксимации.

Однако следует сказать, что применение этих формул в расчетах ограничивает выбор сеточного числа Фурье (половины от числа Куранта). Забегая вперёд, скажем, что приемлемая точность численного решения будет достигнута при для формулы:

. (2.47)

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

, (2.48)

то после решения получим формулу для вычисления значений функции вида:

. (2.49)

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

2.4 Метод неопределённых коэффициентов

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

Рассмотрим теоретическую основу использования данного метода при повышении порядка аппроксимации в МКР и разберём его работу на практическом примере. Для построения формул численного дифференцирования и интегрирования можно использовать метод неопределённых коэффициентов. Заменим производную и значение некоторыми выражениями:

(2.50)

. (2.51)

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

. 2.52)

Уравнению (2.50) соответствует конечно-разностная схема:

. (2.53)

Запишем погрешность аппроксимации исходного дифференциального уравнения схемой:

. (2.54)

Запишем определение. Разностная схема аппроксимирует дифференциальную на отрезке , если выполняется условие:

при . (2.55)

Исходя из того, что и , получим:

. (2.56)

Предположим, что все производные решения до порядка q ограничены:

при (2.57)

Представим все значения и с помощью формулы Тейлора следующим образом:

(2.58)

(2.59)

где ;

.

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

, (2.60)

где ;

;

.

Исходя из уравнения (2,52), имеем:

, (2.61)

где .

Стоит отметить, что производя более аккуратную оценку, можно уменьшить значение в оценке . Если , то и схема (2.54) имеет m-й порядок аппроксимации. Принято считать, что любая разностная схема m-го порядка аппроксимации является схемой q-го порядка аппроксимации при q < m. Если , а , то порядок аппроксимации строго равен m.

Согласно (2.50) и (2.51) для любой гладкой функции справедливы соотношения:

, (2.62)

(2.63)

Запишем лемму для доказательства справедливости уравнений (2.62) и (2.63) в отношении гладких функций: соотношения (2.62) и (2.63) выполнены тогда и только тогда, когда выполняется условие:

. (2.64)

Доказательство. Согласно формуле Тейлора, имеем:

(2.65)

(2.66)

Подставляя соотношения (2.64) для левых частей в формулах (2.62) и (2.63), получаем:

(2.67)

(2.68)

Для справедливости этих соотношений необходимо и достаточно выполнения условий:

(2.69)

(2.70)

(2.71)

Левая часть равенства (2.69) равна, разность левых частей равенств (2.70) и (2.71) равна . Отсюда и следует справедливость утверждения леммы.

Уравнения образуют однородную систему линейных алгебраических уравнений относительно 2k+2 неизвестных. Если число неизвестных больше числа уравнений, т.е. 2k+2 > m+1 (или 2k > m), то эта система имеет ненулевое решение. Можно показать, что при 2k=m эта система имеет однопараметрическую группу ненулевых решений:

(2.72)

. (2.73)

Причем . Выбирая , получим разностную схему 2k-го порядка аппроксимации. Иногда возникает необходимость построить явную схему . Решая систему уравнений , и выбирая решение с , получаем схему (2k-1)-го порядка аппроксимации.

Рассмотрим примеры практического использования вышесказанного. При k = 2 общее решение системы уравнений имеет вид:, , , , , .

Из условия (2.64) (лемма) получаем , а расчётная схема имеет вид:

. (2.74)

Заметим, что если равенства

, (2.75)

(2.76)

выполняются с точностью до членов порядка , то величина , равная разности погрешностей этих соотношений, имеет порядок и, следовательно, согласно (2.53) . Однако для выполнения условия не обязательно требовать, чтобы погрешности соотношений (2.75) и (2.76) имели порядок . К примеру, для последней схемы , в то время как погрешности выражение (2.75) имеет порядок аппроксимации .

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

Пример метода данной группы:

(2.77)

(2.78)

(2.79)

Здесь . Погрешность на шаге у данного метода вычислений .

3. ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ В ЧИСЛЕННЫХ МЕТОДАХ

3.1 Знакомство с устойчивостью в МКР

Первым этапом при обучении студентов любой ступени подготовки анализу устойчивости в численных методах предлагается упрощенное знакомство без глубоких математических доказательств с данной характеристикой.

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

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

...