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

(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) особенно привлекательны при записи в форме трехдиагональной матрицы АХ = В.

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

2.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), полученные методом Кранка- Николсона

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

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

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

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

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

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

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

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



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

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

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

Рисунок 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)

with(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.

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

4. http://mi.mathnet.ru/aa123

5. Cooke G., A weakening of the Euclidean property for integral domains and applications to algebraic number theory. I, J. Reine Angew. Math. 282 (1976), 133-156.

6. Журавлев В. Г., Суммы квадратов над о-кольцом Фибоначчи, Зап. науч. семин. ПО-МИ 337 (2006), 165-190.

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

8. Журавлев В. Г., Одномерные разбиения Фибоначчи, Труды 17-й Международной летней школы-семинара по современным проблемам теоретической и математической физики, Казань, 2005, с. 40-55.с

9. Боревич З. И., Шафаревич И. Р., Теория чисел, 3-е изд., Наука, М., 1985.

10. Дробышевича В. И., Дымникова В. П., Ривина Г. С., Задачи по вычислительной математике [Текст]/ В. И. Дробышевича, В. П. Дымникова, Г. С. Ривина. - М.: Наука, 1988. - 478.

11. Конвей, Дж, Слоэн Н. Упаковки шаров, решетки и группы [Текст]/ Дж. Конвей, Н. Слоэн. - М.: Мир, 1990. - 414.

12. Свободная энциклопедия Википедия [Электронная энциклопедия] // Сетевая энциклопедия Wikipedia. 2000. - http://ru.wikipedia.org/wiki/ - Загл. с экрана.

13. Березин И.С., Жидков Н.П. Методы вычислений [Текст]/ И.С. Березин, Н.П. Жидков. - М.: Физматгиз, 1962. - 507.

14. Мэтьюз, Джон, Г., Финк, Куртис, Д. Численные методы. Использование MATLAB, 3-е издание [Текст]/ Г. Джон Мэтьюз, Финк, Куртис. -- М.: Вильяме, 2001. -- 720.

15. Тихонов А.Н., Самарский А.А. Уравнения математической физики [Текст]/ А.Н. Тихонов, А.А. Самарский. - М.: Наука, 1972.

16. Формалев В.Ф., Ревизников Д.Л. Численные методы [Текст]/ В.Ф. Формалев, Д.Л. Ревизников - М.: ФИЗМАТЛИТ, 2004. - 400.

17. Пирумов У.Г., Численные методы [Текст]/ У.Г. Пирумов. - М.: МАИ, 1998. - 346.

18. Калиткин Н.Н., Численные методы [Текст]/ Н.Н. Калиткин. - М.: Наука, 1976. - 537.

19. Hunt, Brian R Matlab R2007 с нуля [Текст]/ -М.:Лучшие книги, 2008. - 352.

20. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы [Текст]/ Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. - М.:БИНОМ, 2007. - 636.

21. Гречко Л.Г., Сугаков В.И., Томасевич О.Ф., Федорченко А.М. Сборник задач по теоретической физике [Текст]/ Л.Г. Гречко, В.И. Сугаков, О.Ф. Томасевич , А.М. Федорченко. - М.: Высшая школа, 1972. - 336.

22. Демидович Б.П., Марон И.А. Основы вычислительной математики [Текст]/ Б.П. Демидович, И.А. Марон. - М.: Наука, 1966. - 664.

23. Кочин Н.Е. Векторное исчисление и начала тензорного исчисления [Текст]/ Н.Е. Кочин. - М.:Наука,1965. - 425.

24. Мэтьюз Дж., Финк Г., Куртис Д. Численные методы. Использование Matlab [Текст]/ Дж. Мэтьюз, Г. Финк, Д. Куртис. - ИД.:Вильямс, 2001. - 720.

25. Победря Б.Е. Численные методы в теории упругости и пластичности [Текст]/ Б.Е. Победря.- М.: МГУ, 1995. - 366.

26. Поршнев С.В. Вычислительная математика [Текст]/ С.В. Поршнев. - СПб.:БХВ-Петербург,2004. - 320.

27. Срочко В.А. Численные методы [Текст]/ В.А. Срочко. - С.-П.:Лань, 2010. - 208.

28. Бернштейн С.Н. О свойствах однородных функциональных классов // Доклады Ак. Наук СССР,-1947.-№57.-с.111-114.

29. Дзядык В.К. Введение в теорию равномерного приближения функций полиномами. -М.: Наука.-1977.-с.512.

30. Стечкин С.Б. О порядке наилучших приближений непрерывных функций // Изв. АН СССР-Математика,-1931.-№15.-с.219-242.

ПРИЛОЖЕНИЯ

Приложение А

Листинг программы для расчета по методу Кранка-Николсона

f.m

function F=f(x)

F=sin(pi*x)+sin(3*pi*x);

ft.m

function F=f(x,t)

F=sin(pi*x)*exp(-pi^2*t)+sin(3*pi*x)*exp(-9*pi^2*t);

tisys.m

function X=trisys(A,D,C,B)

N=length(B);

for k=2:N

mult=A(k-1)/D(k-1);

D(k)=D(k)-mult*C(k-1);

B(k)=B(k)-mult*B(k-1);

end

X(N)=B(N)/D(N);

for k= N-1:-1:1

X(k)=(B(k)-C(k)*X(k+1))/D(k);

end

crnich.m

function [U,Y]=crnich(c1,c2,a,b,c,n,m)

clc;

% - c1=u(0,t) и c2=u(a,t)

% - а и b - правые точки интервалов [0,а] и [0,Ь]

% - с - постоянная уравнения теплопроводности

% - n и m - число точек решетки на интервалах [0,а] и [0,Ь]

%Выход - U - матрица решений

%Инициализация параметров и матрицы U

h=a/(n-1);

k=b/(m-1);

r=c^2*k/h^2;

s1=2+2/r;

s2=2/r-2;

U=zeros(n,m);

%Граничные условия

U(1,1:m)=c1;

U(n,1:m)=c2;

%Генерирование первого ряда

U(2:n-1,1)=f(h:h:(n-2)*h)';

%Формирование диагональных и не лежащих на диагонали

%элементов А, вектора постоянных В '

%и решение трехдиагональной системы АХ=В

Vd(1,1:n)=s1*ones(1,n);

Vd(1)=1;

Vd(n)=1;

Va=-ones(1,n-1);

Va(n-1)=0;

Vc=-ones(1,n-1);

Vc(1)=0;

Vb(1)=c1;

Vb(n)=c2;

for j=2:m

for i=2:n-1

Vb(i)=U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1);

end

X=trisys(Va,Vd,Vc,Vb);

U(1:n,j)=X';

end

U=U';

%точное решение и определение погрешности

x=0:h:h*(n-1);

t=0:k:k*(m-1);

for i=1:1:n

for j=1:1:m

Y(i,j)=ft(x(i),t(j));

end

end

Y=Y';

pogr=(abs(Y-U));

max(max(pogr))

surf(U)

xlabel('X');

ylabel('t');

zlabel('U(x,t)');

colorbar;

axis([0 n 0 m 0 max(max(U))]);

Приложение B

> solve(1+x+x^2+x^3+x^4);

> 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);

>

> sin(2*Pi/5):=1/4*2^(1/2)*(5+5^(1/2))^(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});

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