Квазирешетки в прикладных задачах обработки цифровой информации

Квазирешетки в прикладных задачах обработки цифровой информации

Министерство образования и науки Российской Федерации

ГОУ Высшего профессионального образования

Северокавказский государственный технический университет

ДИПЛОМНАЯ РАБОТА

Тема Квазирешетки в прикладных задачах обработки цифровой информации

Студентки Васильцовой Ольги Сергеевны

Специальности 230401 «Прикладная математика»

Ставрополь, 2011

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

. РЕШЕТКИ И ИХ СВОЙСТВА

.1 Основы теории сеточных методов

.2 Квазирешетки и их свойства

. КВАЗИРЕШЕТКИ В ПРИКЛАДНЫХ ЗАДАЧАХ ТЕЧЕНИЯ ЖИДКОСТИ

.1Уравнение теплопроводности

.2 Устойчивость. Исследование устойчивости методом гармонического анализа

.3 Постановка задачи

.4 Квазирешетки с применением полиномов Бернштейна

ЗАКЛЮЧЕНИЕ

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

ПРИЛОЖЕНИЯ

ПРИЛОЖЕНИЕ

ВВЕДЕНИЕ

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

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

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

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

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

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

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

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

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

В ходе разработки программы защиты необходимо решить следующие задачи:

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

. Минимизировать время разработки.

. Произвести тестирование разработанной программы.

. Проверить отказоустойчивость программы защиты.

. Минимизировать стоимость разработки.

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

В заключении подведены итоги разработки и проведена оценка успешности разработанного алгоритма.

В приложении содержится листинг расчетных программ.

1. РЕШЕТКИ И ИХ СВОЙСТВА

1.1     Основы теории сеточных методов

Возьмем ограниченную область D из Rⁿ с границей ∂D. Где Rⁿ евклидово пространство n-мерных вещественных векторов. Разобьем пространство Rⁿ на элементарные ячейки {x≡(x₁, x₂,…,x_n): k_ih_i<x_i<(k_i+1)h_i, i=1, 2, …, n}, где k_i - целые числа, h_i=const>0. Величина h_i по переменной x_i, называется шагом сетки. Если h_i=const, то сетка по x_i называется равномерной. Величина h_i называется шагом сетки. Пусть  - объединение элементарных ячеек (включая их границы). Положим, что , а  граница области . Тогда множество вершин ячеек принадлежащих  так же обозначим символом  и назовем его сеткой области , а вершины - узлами сетки. Внутренними узлами сетки,  будем называть вершины ячеек принадлежащие области D. В таком случае граничные узлы составляют множество  и обозначаются .

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

На рисунке 1. и 2 представлены некоторые варианты сеток

Рисунок 1. 1 - Функция одной переменной Ф, заданная на структурированной сетке {x_k}

Рисунок 1. 2 - Функция двух переменных Ф, заданная на структурированной сетке {x_k}

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

Качество построения расчетной сетки в значительной степени определяет успех (неудачу) численного решения уравнения.

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

−       отображение должно быть однозначным;

−       узлы сетки должны сгущаться в области предполагаемой появление больших градиентов;

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

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

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

,

,

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

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

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

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

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

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

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

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

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

 в D,

 на ∂D, (1.1)

где А - линейный оператор, , . В данном случае Ф и F - гильбертовы пространства с областями определения элементов в  и  соответственно, а - линейный оператор граничного условия, , где G - гильбертово пространство функций с областью определения .

Наряду с уравнением (1.1) рассмотрим уравнение в конечномерном пространстве сеточных функций

 в ,

на , (1.2)

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

Введем в сеточных пространствах , G, ,  соответственно нормы , . Пусть  - линейный оператор, который элементу  ставит в соответствие элемент  так, что . При таких условиях, задача (1.2) аппроксимирует задачу (1.1) с порядком n на решении φ, если существуют положительные константы ,  такие, что для всех  выполняются неравенства

(1.3)

.

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

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

 (1.4)

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

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

Рассмотрим задачу:

 в D,  на ∂D (1.5)

Положим, что областью определения D является множество {(х, у): 0 < x < 1, 0 < y < 1}, а f - гладкая функция.

Квадрат  покроем равномерной по х и по у сеткой с шагом h. Узлы области будем отмечать двумя индексами (k, l), где первый индекс k (0 ≤ k ≤ n) соответствует точкам деления по координате x, а второй индекс l (0 ≤ l ≤ n) - по у. Рассмотрим следующие аппроксимации:

где , , ,  - разностные операторы, определенные как

,

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

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

 в ,

 на , (1.6)

где  - множество узлов, принадлежащих границе. С учетом изложенного задача (1.6) может быть приведена к виду

 в ,

 на , (1.7)

где и  - векторы с компонентами и  и

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

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

где - произвольная точка области

Аналогичное разложение будем иметь и для функции f(x, y). Введем в пространстве F_h норму

.

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

, (1.8)

где

.

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

Из последнего утверждения и оценки (1.8) следует, что задача (1.7) аппроксимирует задачу (1.5) со вторым порядком на решениях задачи (1.5), имеющих ограниченные четвертые производные.

.2 Квазирешетки и их своиства

Одномерная квазирешетка Фибоначчи F представляет собою подмножество из R, состоящее из точек  с параметрами N из множества целых чисел Z, где [*] обозначает целую часть числа и  - золотое сечение. Квазирешетка F не является периодической, т.е. равенство F +t = F возможно только при t = 0.Если поставить себе целью найти другие аналогичные F квазирешетки, то один из естественных подходов - это ввести в рассмотрение квадратную квазирешетку Фибоначчи F² = F × F ( рис. 1.1). С помощью „Cut and Project" метода квазирешетку F² можно характеризовать равенством

(F²)' = O² ∩ J², (1.9)

где О - кольцо целых чисел квадратичного поля , ' - сопряжение в F и  [2]. Здесь в равенстве (1.9) J² = J × J - «окно» или внутреннее пространство, выделяющее квазирешетку F² из всюду плотного множества O² в виде (r, R) - системы Делоне, т.е. не слишком плотной и не слишком разреженной системы точек на плоскости R².

Рисунок 1.1 - Квадратная решетка Фибоначчи F²

Двумерные квазирешетки F² первым ввел и исследовал их физические свойства Лифшитц [3]. В данном случае F² используеться как двумерная квазирешетка, которая содержит в себе все основные виды одномерных квазирешеток Фибоначчи, что и позволяет получить их полную классификацию [4].

Выберем две произвольные точки из F² и проведем через них прямую L. Тогда данная прямая L порождает

(1.10)

одномерную квазирешетку L, вложенную в F². В настоящей работе классифицируются все возникающие указанным способом квазирешетки L из F². В основу положена параметризация

(1.11)

квазирешетки L точками t' из множества параметров , где γ - некоторый интервал из J,  - примитивный направляющий вектор прямой L из уравнения (1.10) и υ - вектор сдвига.

Одномерные квазирешетки  и  называются подобными , если одна переводится в другую некоторым преобразованием подобия плоскости R². Т. е. две квазирешетки L, L ' из квадратной квазирешетки F² подобны:

. (1.12)

Вторая эквивалентность из (1.12) означает , где  - подгруппа положительных единиц из группы единиц О^× кольца Фибоначчи О. Используя параметризацию (1.11) можно показать, что множество классов подобных квазирешеток находится во взаимно-однозначном соответствии с числами из множества .

Две одномерные квазирешетки L, L ' из F² принадлежат к одному локальному типу или локально эквивалентны , если найдется такое преобразование подобия s, что квазирешетки  имеют одно и то же множество расстояний между соседними точками. Это множество может состоять из двух g, e или трех g, e, g + e расстояний, где g/e = τ. Выберем произвольную точку x пусть x_-, x₊ ϵ L - соседние для x точки, упорядоченные, согласно направляющего вектора β. Тогда найдется такая точка x' ϵ l', для которой соседние расстояния будут , и наоборот. Из определения следует, что локально эквивалентные квазирешетки l, l^' локально не различимы.

В таблице. 1.1 для квазирешеток l локального типа j = 0,1, 2, 3 указаны все допустимые первые локальные окружения точек. Пусть, например, j = 0 и точка x ? L имеет соседние точки x_- ,x₊ ? L. Тогда расстояния  могут принимать следующие три значения: (r_-, r₊) = (g, e), (e, g) или (e ,e). В таблице 1.1 можно видеть, как при изменении длины интервала γ для квазирешетки  появляются или исчезают те или иные виды локальных окружений. Во втором столбце указано количество всех локальных окружений. Локальные типы j =0, 1, 2, 3 различны, а следующий локальный тип j =0^' эквивалентен j =0. При переходе с одного уровня на другой происходит увеличение локальных расстояний в γ + 1 раз, что влечет периодичность j mod 4.

Таблица 1.1 - Локальные типы одномерных квазирешеток L из F²

Также известно, что любая одномерная квазирешетка l из двумерной квазирешетки Фибоначчи F² имеет один из локальных типов j, где j =0, 1, 2, 3, и каждый такой локальный тип реализуется некоторой квазирешеткой l из F².

Любые две решетки одной размерности над Z подобны. Данный факт резко контрастирует с существованием счетного множества классов подобных одномерных квазирешеток l из F². Более того, даже на локальном уровне существуют четыре различных типа таких квазирешеток. Все они образуют класс одномерных квазирешеток Фибоначчи, которые характеризуются как квазирешетки, параметризуемые вращениями окружности и их индуцированными отображениями (отображениями Пуанкаре).

Уровень m = 0,1, 2,... квазирешетки  определяется из соотношений

 или

и соответственно квазирешетка L называется невырожденной или вырожденной. В первом случае интервал γ разбивается

(1.13)

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

Разобьем  множество точек квазирешетки L имно-жество ее параметров  аналогично разбиению (1,13). Предположим, что β ≠ 0 для направляющего вектора β ≠ (β₁, β₂ ),и пусть p r₁x = x₁ - проекция вектора x = (x₁, x₂) на ось OX .Для каждого вида точек x = x_Y из L , где Y = G, GE, E, определим частоту

его появления на квазирешетке . Если квазирешетка имеет четный уровень m, то частоты py ее точек x_Y ϵ L вычисляются по формулам

(1.14)

если L - невырожденная квазирешетка, и

 (1.15)

если L вырожденная. Аналогичные формулы имеют место для нечетного уровня m. Пусть n_L(X) обозначает количество точек x = (x₁, x₂) на квазирешетке L с 0 ≤ x₁ ≤ X. Тогда для n_L(X) имеет место асимптотическое равенство

 при (1.16)

где .

Кольцо Фибоначчи  является евклидовым. Оно обладает алгоритмом деления с остатком по норме [1]. Еще одно применение одномерных квазирешеток состоит в том, что нормированная квазирешетка  задает, по терминологии Кука, одномерный 2-ступенчатый алгоритм Евклида, когда в качестве неполного частного для , где  и , выбирается максимальное число вида

с номером N = 0,1, 2,... Получена верхняя оценка для количества шагов в алгоритме, из которой следует, что данный одномерный алгоритм является быстрым[5].

Определим отображение

(1.17)

обладающее свойствами , где , функция r(N) принимает значения из N = {0,1, 2,...}, J = [-1, τ) и черта обозначает замыкание множеств. Тогда  [6]. Отображение можно записать в явном виде

(1.18)

Одномерная квазирешетка Фибоначчи F - это дискретное подмножество из R. Неопределяемое следующим образом (1.18):

,(1.19)

где  для всех N ϵ Z. Заметим, что , при этом ' - сопряжение в квадратичном поле . Квазирешетка F обладает целым рядом интересных свойств. Так, например, успешно применяют ее к построению алгоритма деления с остатком в кольце O . Здесь же приведены некоторые простейшие свойства квазирешетки F.

Обозначим F⁺ подмножество из F, состоящее из точек вида δ(N) с параметрами N ? N. Разностная функция  равна , если , и , если  Отсюда заключаем, что множество F⁺ параметризуется вращениями окружности. Из [7] следует, что последовательность

(1.20)

где , является последовательностью Штурма. Последовательность и (1.20) полностью характеризует положительную часть F⁺ квазирешетки F, так как , если , и , если . Последовательность u инвариантна относительно подстановки, действующей на символы 0, 1 по правилу, и получается из 1 многократным применением подстановки π. Длина последовательностей  равна , где  - числа Фибоначчи, определяемые рекуррентным соотношением F_m+2 = F_m+1 + F_m иначальным условием F₁ = 1, F₂ = 2. Так же, как и для чисел Фибоначчи F_m, для последовательностей  выполняется свойство рекуррентности

(1.21)

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

(1.22)

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

Множество  обладает центральной симметрией  для любого N ϵ Z\ {-1}Используя формулу (1.18), можем записать i на языке параметров . Отсюда получаем, что проколотая квазирешетка  центрально-симметрична

с центром симметрии , не являющимся элементом из .

Кроме симметрии i, квазирешетка F допускает богатую полугруппу подобий ,  - подмножество тех N ? Z, для которых |. Для любого N ?  выполняется включение

.

Особый интерес представляет случай N =1, когда наблюдается подобие

(1.23)

с коэффициентом , являющимся основной единицей кольца Фибоначчи . В данном случае подобие (1.23) представляет собою инфляцию проколотой квазирешетки , соответствующую подстановке π.

Квазирешетка F не является периодической, т.е. равенство F + t = F возможно только при t = 0. Это вытекает из аналогичного свойства множества параметров . Однако имеют место следующие два свойства:

 (1.24)

где , означающие, что квазирешетка F «почти» замкнута относительно операций сложения и вычитания (свойство Мейера).

С целью расширить класс одномерных квазирешеток введем квадратную квазирешетку (рисунок 1.1)

(1.25)

Данную квазирешетку можно характеризовать как множество всех точек (x₁, x₂) из O² с ограничением (x’₁, x’₂) ϵ J² или более кратко -

.(1.26)

Пусть прямая L проходит через две точки из F². Тогда будет одномерной квазирешеткой. В частности, если L - главная диагональ, то получаем квазирешетку L, подобную квазирешетке F .А если L - дополнительная диагональ, то L уже не будет подобна F, поскольку расстояния между соседними точками в L принимают три различных значения. Как например у квазирешетки L_2,0 на рисунке 1.1.

Чтобы классифицировать возникающие таким способом одномерные квазирешетки L, нам потребуется детальное изучение разбиений Фибоначчи полуинтервала J.

По существу мы используем известный метод «Cut and Project», выбирая в определении (1.26) в качестве «окна» или внутреннего пространства W квадрат J². Если выбрать W произвольным выпуклым множеством, то новых типов квазирешеток L не добавится [2]. Если же допустить невыпуклые W, то могут возникнуть квазирешетки L получающиеся объединением предыдущих квазирешеток.

Отождествим полуинтервал I =[0,1) с единичной окружностью R/Z, и определим нанем сдвиг

,

где  - золотое сечение. Поскольку сдвиг τ и 1 несоизмеримы, то можно для любого ε ϵ I определить индуцированное отображение (или отображение Пуанкаре)

, где  (1.27)

интервал, на котором действует S_ε. Известно, что индуцированное отображение S_ε. изоморфно 1) нециклическому перекладыванию трех полуинтервалов из  или 2) перекладыванию двух полуинтервалов из , что эквивалентно вращению единичной окружности на некоторый угол. Случай 2) возникает, когда параметр ε принимает значения

 для (1.28)

Для таких е полуинтервал  разбивается

(1.29)

на два полуинтервала  с соотношением длин  и  для m четных и нечетных соответственно. В (1.29) обозначает некоммутативную операцию прикладывания одного полуинтервала к другому. Индуцированное отображение  перекладывает

(1.30)

полуинтервалы из разбиения (1.29). При этом индуцированное отображение S^(m) изоморфно

(1.31)

сдвигу S или ему обратному сдвигу S^-1 для четных или нечетных m [8].

Собственное разбиение Фибоначчи  единичного полуинтервала I уровня m получается, если подействовать несколько раз сдвигом S на полуинтервалы  до их первого возвращения в полуинтервал 1т.Разбиение  содержит F_m+1 и F_m+2 полуинтервалов длины и  и ,если уровень m четный, и соответственно F_m+1 и F_m+2 полуинтервалов, если m нечетный. Тогда инвариантность разбиений  относительно действия сдвига S:

(1.32)

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

.3 Алгоритмы построения квазирешетки и оценка их свойств

При переходе от неоднородной квазирешетки  к отвечающей ей однородной квазирешетке  необходимо знать вектор сдвига υ. Для этого требуется уметь находить частное решение (x₁₀,x₂₀) из O² неоднородного уравнения

(1.32)

Кольцо Фибоначчи O является евклидовым относительно деления с остатком по норме, следовательно, частные решения уравнения (1.32) можно найти с помощью алгоритма Евклида. Отметим, что указанный алгоритм весьма трудоемкий [9].

С другой стороны, использование одномерной квазирешетки F (1, 3) позволяет в кольце Фибоначчи O определить деление с остатком полностью аналогично тому, как это делается в кольце целых рациональных чисел Z, и тем самым получить быстрый алгоритм деления в кольце O. Возникающий таким образом алгоритм Евклида является 2-ступенчатым и, что необходимо отметить, по своей сути является одномерным. В теореме 1.1 доказана верхняя оценка для количества шагов в алгоритме [5].

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

 где .

Необходимость замены основной квазирешетки F сжатой квазирешеткой F₁ обусловлена тем, что разностная функция

 или τ (1.34)

принимает положительные значения, не превосходящие 1 Данное обстоятельство позволяет в алгоритме деления с остатком формально использовать квазирешетку F₁ вместо кольца целых рациональных чисел Z.

Таблица 1.2 - Замены основной квазирешетки F сжатой квазирешеткой F₁

Пусть , где . Из (1.34) следует существование единственного , удовлетворяющего условиям

, где (1.35)

Используя (1.35), получаем представление

, где  и .(1.36)

₁-алгоритм - это рекуррентное повторение действия (1.36):

(1.37)

где  для i = 0,1, 2,...

Докажем конечность F₁ - алгоритма (1.37), т.е. покажем, что  для некоторого . Пусть ,тогда имеем

.(1.38)

Если N_i ≥ 4, то получаем неравенства

,(1.39)

так как из включения  следует .

Аналогично имеем

 для ;

 для ;(1.40)

 для .

Итак, в (1.39) для  мы видим препятствие, не позволяющее доказать конечность F₁ - алгоритма. Это означает, что (1.37) не является обычным алгоритмом Евклида. И действительно, далее мы убедимся в том, что (1.37) - это 2-ступенчатый алгоритм по терминологии Кука [5]. Для доказательства потребуется дополнительная информации нам потребуется дополнительная информация о величине  на следующем шаге .

Если , то ввиду (1.40) получаем неравенство

,(1.41)

так как , иначе в силу (1.38) выполнялось бы неравенство

 или . Откуда следует , что противоречит предположению . Поэтому для имеем (таблица 1.2)

.(1.42)

Пусть , тогда, снова используя (1.40), приходим к неравенству

,

из которого и (1.40) для  получаем

.(1.43)

Аналогично в остальных случаях  имеем

и, следовательно,

.(1.44)

Поскольку , то можем переписать F₁ - алгоритм (1.37) в более привычном виде

(1.45)

Оценим сверху величину . Применяя несколько раз равенство , получаем для любого i ≥0 связь

(1.46)

между остатками r_i+i и r₀. Удобно указанную связь (1.46) перевести в матричную форму. Так как , то отсюда приходим к следующему матричному соотношению

,(1.47)

где M_i = Q_o…Q_i и матрицы Q_j имеют вид . Пусть  обозначает длину вектора . Тогда, очевидно,

(1.48)

и правая часть, ввиду соотношения (1.47), может быть оценена как

(1.49)

.(1.50)

Пусть υ - вектор из R²,  с коэффициентом r ≥ 0, (•, •) - скалярное произведение в R². Воспользуемся неравенством

,(1.51)

вытекающим из следующей цепочки преобразований , где Λ - максимальное собственное значение матрицы RR^t. Поскольку R = R^t - симметрическая матрица, то

, где (1.52)

максимальное собственное значение матрицы R. Теперь из (1.51), (1.52) получаем нужное нам неравенство

.

Применим его к правой части (1.49). Имеем

(1.53)

(1.54)

собственное значение матрицы R_j (1.50). Отсюда и (1.48) получаем неравенство , из которого и (1.46) вытекает неравенство для нормы остатков r_i+1 в .F₁-алгоритме (1.45)

.(1.55)

Согласно (1.39), (1.40) и (1.54), получаем следующие неравенства для произведений θ_i λ_i:

 для

 для (1.56)

 для .

Отсюда, используя (1.39) и (1.40), оцениваем сверху попарные произведения , предварительно перегруппировывая их к виду :

 для , ;

 для , ;

 для , .(1,57)

Разобьем произведение  в правой части неравенства (1.55) на 1) сдвоенные пары  для ; и 2) одиночные пары  для N_k ≥ 2 и не входящих в 1). Оценим их с помощью неравенств (1.56), (1.57) и применим к (1.55). Получаем

, если i четное,(1.58)

 если i нечетное.

Теперь отметим, что в F₁-алгоритме (1.45) остатки r_i+1 принадлежат кольцу O,следовательно, , если r_i+1 =0. Отсюда и из неравенств (1.58) вытекает

Теорема 1.1. Для любых x, y из кольца Фибоначчи , F₁-алгоришм (1.45) заканчивается не более чем через

12(ln(x) + ln(|(x', y')|)) + 1(1.59) шагов.

Замечание Численные вычисления показывают, что множитель 12 в (1.59), как правило, можно заменить на 1.

2. КВАЗИРЕШЕТКИ В ПРИКЛАДНЫХ ЗАДАЧАХ ТЕЧЕНИЯ ЖИДКОСТИ

математика квазирешетка полином вычислительный

2.1 Уравнение теплопроводности

Допустим, что рассматривается некоторое тело и изучается его тепловое состояние. Последнее будет известно, если для каждой точки тела мы будем знать температуру Т в любой момент; иными словами, тепловое состояние тела характеризуется скалярной функцией T(r,t)=T(x,y,z,t).

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

Рассмотрим внутри тела некоторый объем V, ограниченный поверхностью S, и подсчитаем двумя способами изменение количества тепла, заключенного в объеме V (рисунок 5)

Рисунок 2.1 - Тело с некоторым объемом, ограниченным поверхностью S

Плотность тела обозначим через ρ (если тело неоднородное, то ρ будет функцией точки ρ(r)=ρ(x,y,z)), а теплоемкость через с (в случае неоднородности тела с тоже есть функция точки) [9].

Рассмотрим элемент dV объема; масса этого элемента равна ρdV; за время dt этот элемент нагревается на  градусов; на это требуется, по самому определению теплоемкости, количество тепла, равное

интегрируя по всему объему, увидим, что за время dt всему объему V необходимо было сообщить количество тепла, равное [2].

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

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

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

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

   (2.1)

Если на границах и заданы значения искомой функции в виде

 , , (2.2)

 , , (2.3)

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

 , , (2.4)

то задачу (2.1) - (2.4) называют первой начально-краевой задачей для уравнения теплопроводности (2.1).

В терминах теории теплообмена - распределение температуры в пространственно-временной области

 a² - коэффициент температуропроводности, а (2.2), (2.3) с помощью функций , задают температуру на границах  и .

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

 , , (2.5)

 , , (2.6)

т.е. граничные условия второго рода, то задачу (2.1), (2.5), (2.6), (2.4) называют второй начально-краевой задачей для уравнения теплопроводности (2.1). В терминах теории теплообмена на границах в этом случае заданы тепловые потоки.

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

 , , (2.7)

 , , (2.8)

т.е. граничные условия третьего рода, то задачу (2.1), (2.7), (2.8), (2.4) называют третьей начально-краевой задачей для уравнения теплопроводности (2.1). В терминах теплообмена граничные условия (2.7), (2.8) задают теплообмен между газообразной или жидкой средой с известными температурами на границе и  на границе  и границами расчетной области с неизвестными температурами , . Коэффициенты α, β - известные коэффициенты теплообмена между газообразной или жидкой средой и соответствующей границей.

Для пространственных задач теплопроводности в области  первая начально-краевая задача имеет вид

 (2.9)

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

Согласно методу сеток в плоской области D строится сеточная область D_h, состоящая из одинаковых ячеек. При этом область D_h должна как можно лучше приближать область D. Сеточная область (то есть сетка) D_h состоит из изолированных точек, которые называются узлами сетки. Число узлов будет характеризоваться основными размерами сетки h: чем меньше h, тем больше узлов содержит сетка. Узел сетки называется внутренним, если он принадлежит области D, а все соседние узлы принадлежат сетке D_h. В противном случае он называется граничным. Совокупность граничных узлов образует границу сеточной области Г_h.

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

Нанесем на пространственно-временную область ,  конечно разностную сетку ω_h,τ:

 (2.10)

с пространственным шагом h=l/N и шагом по времени τ=T/K.

Рисунок 2.2 - Конечно-разностная сетка

Введем два временных слоя: нижний ,на котором распределение искомой функции u(x_j,t^k), , известно (при к = 0 распределение определяется начальным условием (2.4) u(x_j,t^k)=ψ(x_j)), и верхний временной слой t^k+1=(k+1) τ, на котором распределение искомой функции u(x_j,t^k+1), .

Сеточной функцией задачи (2.1) - (2.4) называют однозначное отображение целых аргументов j,k в значения функции .

На введенной сетке вводят сеточные функции ,  первая из которых известна, вторая подлежит определению. Для определения в задаче (2.1) - (2.4) заменяют (аппроксимируют) дифференциальные операторы отношением конечных разностей (более подробно это рассматривают в разделах численных методов «Численное дифференцирование»), получают

, (2.11)

, (2.12)

Подставляя (2.11), (2.12) в задачу (2.1) - (2.4), получим явную конечно-разностную схему для этой задачи в форме

 (2.13)

В каждом уравнении этой задачи все значения сеточной функции известны, за исключением одного, , которое может быть определено явно из соотношений (2.13). В соотношения (2.13) краевые условия входят при значениях j=1 и j=N-l, a начальное условие - при k = 0.

Если в (2.12) дифференциальный оператор по пространственной переменной аппроксимировать отношением конечных разностей на верхнем временном слое:

, (2.14)

то после подстановки (2.11), (2.14) в задачу (2.1) - (2.4), получим неявную конечно-разностную схему для этой задачи:

 (2.15)

Теперь сеточную функцию на верхнем временном слое можно получить из решения (15) с трехдиагональной матрицей. Эта СЛАУ в форме, пригодной для использования метода прогонки, имеет вид

;

 ;

;

;

;

;

.

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

Рисунок 2.3 - Шаблон явной конечно-разностной схемы для уравнения теплопроводности

Рисунок 2.4 - Шаблон неявной конечно-разностной схемы для уравнения теплопроводности

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

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

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

2.2 Устойчивость. Исследование устойчивости методом гармонического анализа

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

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

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

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

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

, (2.16)

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

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

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

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

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

 (2.17)

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

 (2.18)

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

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

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

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

 (2.19)

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

, (2.20)

т. е. условие (2.20) является необходимым условием устойчивости разность

Явная конечно разностная схема, записанная в форме

 (2.21)

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

 (2.22)

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

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

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

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

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

Рисунок 2.5 - Двусторонний метод аппроксимации

Проведенный анализ дал блестящий пример так называемых двусторонних методов, исследованных В. К. Саульевым

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

 (2.23)

Где θ - вес неявной части конечно-разностной схемы,

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

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

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

,

 (2.24)

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

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

  (2.25)

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

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

Рассмотрим порядок аппроксимации неявно-явной схемы с весами, для чего разложим в ряд Тейлора в окрестности узла (x^j,t_k) на точном решении значения сеточных функций  по переменной t, ,  по переменной х и полученные разложения подставим в (2.23):

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

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

,

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

Используем в уравнение (2.23) подстановку r=a^2k/h². Но в то же время его нужно решить для трех "еще не вычисленных" значений , , и . Это возможно, если все значения перенести в левую часть уравнения. Затем упорядочим члены уравнения (2.23) и в результате получим неявную разностную формулу

 (2.26)

для i=2,3,…,n-1. Члены в правой части формулы (2.26) известны. Таким образом, формула (2.26) имеет вид линейной трехдиагональной системы АХ=В. Шесть точек, используемых в формуле Кранка-Николсона (2.26), вместе с промежуточной точкой решетки, на которой основаны численные приближения, показаны на рисунке 5.

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

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

, (2.27)

для i=2,3,…,n-1. Граничные условия используются в первом и последнем уравнениях (т. е. в  и  соответственно).

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

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

.3 Постановка задачи

Используем метод Кранка-Николсона, чтобы решить уравнение

,

где x?(0;1),

t?(0;0.1),

с начальным условием

,

где t=0,

x?(0;1),

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

u(0,t) = g₁(t) ≡ 0,

u(1,t) = g₂(t) ≡ 0.

Решение будем искать в ППП MatLab 2007. Создадим четыре выполняемых m-фала: crnich.m - файл-функция с реализацией метода Кранка-Николсона; trisys.m - файл-функция метода прогонки; f.m - файл-функция задающая начальное условие задачи; fе.m - файл-функция задающая функцию определяющую точное решение задачи(найдена аналитическим путем). Листинги программ представлены в приложении А.

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

Решения для данных параметров отразим в таблице 1. Трехмерное изображение данных из таблицы покажем на рисунке 5.

Таблица 2.1 - Значения u(х_i, t_i), полученные методом Кранка- Николсона

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

u(x,t) = sin(πx)e^-π2t+ sin(3πx)e^-9π2t,

истинные значения для последнего представлены в таблице 2.2

Максимальная погрешность для данных параметров равна 0,005

Таблица 2.2 - точные значения u(х_i, t_i), при t=0.1

Рисунок 2.7 - Решение u=u(х_i, t_i), для метода Кранка-Николсона

.4 Квазирешетки с применением полиномов Бернштейна

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

Базисные многочлены Бернштейна степени n находятся по формуле

и обладают следующими свойствами:

 при  и всех k=0, 1,…,n

Рисунок 2.8 - Заданная область

Опишем, заданную область уравнениями границы:

.

И нормали к данным уравнениям выглядят так:

.

Тогда уравнения границ будут выглядеть следующим образом:

Запишем коэффициенты обобщения полинома Берштейна с учетом того, что

:

При n=1 узлы решетки будут находиться в вершинах пятиугольника. При n=2, к узлам в вершинах прибавляется еще узлы на серединах сторон и внутренний узел в данном случае это будет точка с координатой (0.3808, 0.5401) (приложение B).

Для нахождения узлов решетки и значения функции в этой точке, для случая n=3, решаются дифференциальные уравнения вида:

с начальными условиями

Применяя в Maple 11 такую последовательность операторов:

Рисунок 2.9 - Нахождение внутреннего узла заданной области

fsolve({diff(ld1^2*ld3,x)=0,diff(ld1^2*ld3,y)=0},{x=0.2, y=0.5});

evalf(ld1*ld2);

(1. + 0.8090169942 x - 0.8312538757 y) (1. + 1.144122805 x + 0.8312538757 y)^2(1. - 0.3090169942 x + 0.5877852520 y)^2*U^2*(1. - 1. x)(DEtools):

> DE := (diff(U, x))^2+(diff(U, y))^2;

> with(Optimization);

> Minimize(DE);

ЗАКЛЮЧЕНИЕ

Основным результатом дипломной работы является разработка алгоритма разбиения квазирешеткой и его программный расчет. В качестве среды разработки программного средства была выбрана MatLab R2007b и Maple 11 в связи с тем, что они позволяет качественно и быстро произвести расчет.

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

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

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

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

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

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

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

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

1.       Марчук, Г.И. Методы вычислительной [Текст]/ Г.И. Марчук. - М.: Наука, 1989. - 608.

2.       Moody R. V., Model sets: a survey, Quasicrystals to More Complex Systems (F. Alex, F. Denoyer, and J. P. Gazeau, eds.), EPD Science, Les Ulis, and Springer-Verlag, Berlin, 2000, pp. 145-166.

.        Lifshitz R., The square Fibonacci tiling, J. Alloys Compounds 342 (2002), 186-190.

4.       <#"539241.files/image432.gif">

> evalf(cos (2*Pi/5));

> evalf([-1/4+1/4*5^(1/2)+1/4*I*2^(1/2)*(5+5^(1/2))^(1/2), -1/4-1/4*5^(1/2)+1/4*I*2^(1/2)*(5-5^(1/2))^(1/2), -1/4-1/4*5^(1/2)-1/4*I*2^(1/2)*(5-5^(1/2))^(1/2), -1/4+1/4*5^(1/2)-1/4*I*2^(1/2)*(5+5^(1/2))^(1/2)]);

> 

> cos(2*Pi/5):=-1/4+1/4*5^(1/2);

> 

> cos(4*Pi/5):=-1/4-1/4*5^(1/2);

> sin(Pi/5):=(1-cos(2*Pi/5))^(1/2);

> sin(2*Pi/5):=1/4*2^(1/2)*(5-5^(1/2))^(1/2);

> cos(Pi/5):=(1+cos(2*Pi/5))^(1/2);

> W0:=1-x;

> W1:=1-cos(2*Pi/5)*x-sin(2*Pi/5)*y;

> W2:=1-cos(4*Pi/5)*x-sin(4*Pi/5)*y;

> W3:=1-cos(6*Pi/5)*x-sin(6*Pi/5)*y;

> W4:=1-cos(8*Pi/5)*x-sin(8*Pi/5)*y;

> W3 := 1+cos(1/5*Pi)*x+sin(1/5*Pi)*y;

> W2 := 1-(-1/4-1/4*5^(1/2))*x-sin(1/5*Pi)*y;

> W4 := 1-cos(2/5*Pi)*x+sin(2/5*Pi)*y;

> ld1:=W2*W3*W4*U;

> ld2:=W0*W3*W4*U;

> ld3:=W0*W1*W4*U;

> ld4:=W2*W0*W1*U;

> ld5:=W2*W3*W1*U;

> U:=1/(ld1+ld2+ld3+ld4+ld5);

> simplify(1/U);

> plot3d(ld1^2*ld3,x=-1..1, y=-1..1,view=0..0.2);

> fsolve({diff(ld1^2*ld3,x)=0,diff(ld1^2*ld3,y)=0},{x=0.2, y=0.5});