Математические модели технических объектов

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ ПОЛИТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

ФАКУЛЬТЕТ ТЕХНИЧЕСКОЙ  КИБЕРНЕТИКИ

КАФЕДРА "ИНФОРМАЦИОННЫХ И УПРАВЛЯЮЩИХ СИСТЕМ"

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Отчет

по курсовой работе

по предмету “Математические модели

технических объектов”

 

 

 

 

 

 

 

 

 

Выполнила:

Васильева Н.В.

Гр. 3084/3

Проверил:

доц. С.П.Воскобойников.

 

 

Содержание

 

 

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

 

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

 

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

 

          

 

Вид области  и расположение граничных условий: 

 

 

 

 

 

 

 

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

 

 

 

 

 

 

 

 

 

 

 

 

Дискретная модель.

 

Построим дискретную модель для данной двумерной модели.

 

a=x_L<x₁< ...<x_R=b

c=y_L<y₁< ... <y_R=d

 

Введем в прямоугольнике [x ,x ]´[y ,y ] равномерную для простоты сетку, то есть множество узлов с координатами

                                                          (2)

где ,  N и N - числа разбиений по x и по y. Кроме этой основной сетки, введем вспомогательную сетку с координатами .

   

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

Интегрируем исходное ДУ второго порядка в промежутке , и, аппроксимируя полученные производные,  получаем систему алгебраических уравнений для всех внутренних точек области:

 

для i=1,2,..., N -1;  j=1,2,..., N -1:

 

 (1)

 

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

Для левой вертикальной границы (х=x_L).

i=0, j=1..N_y-1

(2)

 

Для границы y=y_L:

i=1,2,..., N -1;  j=0

  (3)   

Теперь рассмотрим вид  алгебраического уравнения (1) для  угловых точек области. В этом случае на граничные условия  заменяются две производные.

 

Для x=x0 y=y0:

 

i=0, j=0

 

 (5)

 

Теперь рассмотрим предпоследнее  разбиение по Х (i=N_x-1) и по Y(j=Ny-1). Здесь (i+1)-я и (j+1)-я компоненты в алгебраическом уравнении заменяются на граничные условия первого рода – известные функции V2(y) и V4(x) соответственно.

 

Для решения системы (1) - (5) перенумеруем узлы разбиения по оси Х. Для этого сделаем следующую замену:

 

V_i,j = W_k

V_i-1,j = W_k-1

V_i+1,j = W_k+1

V_i,j-1  = W_k-L

V_i,j+1 = W_k+L,  где L = N_x

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

 

k =  j*N_x + i,

 

где i = 0..N_x-1;   j = 0..N_y-1.

 

При этом размерность системы будет задаваться так:        

N=N_x*N_y.

 

 

 

Решение системы неявным  методом сопряженных градиентов.

 

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

 

A*X=F

следующим образом:

 

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

 

, где  - задаваемая точность решения.

 

Матрица А будет квадратная, пятидиагональная (остальные элементы нулевые), симметричная и имеет следующую структуру:

 

 

      А =

 

 

 

 

 

 

Вследствие такого вида матрицы А будем хранить только три диагонали: А_k, А_k+1 и А_k+L.

 

 

 

При этом в методе сопряженных градиентов:

 

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

N=N_x*N_y

Соответствующие элементы диагоналей Ak, Bk и Ck и b для каждого узла k (i и j):

 

i=1..N_x-1,  j=1..N_y-1

i=0, j=1,..Ny-1

j=0, i=1,..Nx-1

 

 

i=0, j=0

Тесты для заданной модели

Во всех тестах вектор начальных приближений возьмем нулевой.

Тест №1

Первый тест определяет решение интегро - интерполяционного метода без погрешности и имеет вид констант:

 

U(x,y) = 3;

f(x,y) = 0;

g₁(y) = 3;

g₂(y) = 3;

g₃(x) = 3;

g₄(x) = 3;

k₁(x,y) = 2; k₂(x,y) = 5

; ;

 

 

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

k – количество итераций	Nx	Ny	N- порядок матрицы	Погрешность
7	4	4	16	7,46446529475975E-07
12	8	8	64	4,32180643716862E-07
22	16	16	256	9,01495928129492E-07
43	32	32	1024	6,22035156849776E-06
52	64	64	4096	2,78118248218107E-05

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

Тест №2

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

 

U(x,y) = 2*x+3*y+5;

f(x,y) = -10;

g₁(y) = kappa1*(2*xLeft+3*y+5) - (2*xLeft+y+1)*2;

g₂(y) = (2*xRight+3*y+5);

g₃(x) = kappa3*(2*x+3*yLeft+5)-(x+2*yLeft+1)*3;

g₄(x) = 2*x+3*yRight+5;

k₁(x,y) = 2*x + y + 1;

k₂(x,y) = x + 2*y + 1;

; ;

 

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

k – количество итераций	Nx	Ny	N- порядок матрицы	Погрешность
7	4	4	16	1,08685551136745E-06
11	8	8	64	3,01398424440436E-06
17	16	16	256	2,78906810855517E-05
34	32	32	1024	1,84274984151855E-05
64	64	64	4096	3,93771196849357E-05

В данном примере хорошо видно, что функция выходит на свое стационарное значение.

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

Тест №3

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

Зададим:

U(x,y) = 2*x*x + 3*y*y*y + 4;

k₁(x,y) = 1 + 2*x*x + y;

k₂(x,y) = 1 + x + 2*y*y;

; ;

 

Получаем следующие  функции:

f(x,y) = -(24*x*x + 18*x*y + 22*y + 72*y*y*y + 4);

g₁(y) = kappa1*(2*xLeft*xLeft+3*y*y*y+4)-(1+2*xLeft*xLeft+y)*(4*xLeft);

 g₂(y) = 2*xRight*xRight+3*y*y*y+4;

g₃(x) = kappa3*(2*x*x+3*yLeft*yLeft*yLeft+4)-(1+x+2*yLeft*yLeft)*(9*yLeft*yLeft);

g₄(x) = 2*x*x+3*yRight*yRight*yRight+4;

Nx = 16, Ny = 16; N = 256;

Полученные результаты:

 

x	y	U(x,y)	x	y	U(x,y)	x	y	U(x,y)	x	y	U(x,y)
0	0	4,042258511	0,5	0	4,5502058	1	0	6,0413218	1,5	0	8,5248301
0	0,125	4,048213087	0,5	0,125	4,5547821	1	0,125	6,0445089	1,5	0,125	8,527125
0	0,25	4,089864183	0,5	0,25	4,5950917	1	0,25	6,0838567	1,5	0,25	8,5662101
0	0,375	4,201847613	0,5	0,375	4,7059665	1	0,375	6,1942096	1,5	0,375	8,6766218
0	0,5	4,418977791	0,5	0,5	4,9222442	1	0,5	6,4103746	1,5	0,5	8,893061
0	0,625	4,776198033	0,5	0,625	5,278838	1	0,625	6,7672149	1,5	0,625	9,2503374
0	0,75	5,30855187	0,5	0,75	5,8107303	1	0,75	7,2996668	1,5	0,75	9,7834123
0	0,875	6,051153307	0,5	0,875	6,5529811	1	0,875	8,042711	1,5	0,875	10,527309
0	1	7,039152719	0,5	1	7,5407086	1	1	9,0314121	1,5	1	11,517097
0	1,125	8,307727157	0,5	1,125	8,8090593	1	1,125	10,300863	1,5	1,125	12,787866
0	1,25	9,892064806	0,5	1,25	10,393201	1	1,25	11,886176	1,5	1,25	14,374713
0	1,375	11,82735741	0,5	1,375	12,328305	1	1,375	13,8225	1,5	1,375	16,312761
0	1,5	14,14880884	0,5	1,5	14,649561	1	1,5	16,144999	1,5	1,5	18,637124
0	1,625	16,89161685	0,5	1,625	17,392163	1	1,625	18,888828	1,5	1,625	21,382928
0	1,75	20,09096519	0,5	1,75	20,591308	1	1,75	22,089157	1,5	1,75	24,58527
0	1,875	23,78203694	0,5	1,875	24,282193	1	1,875	25,781154	1,5	1,875	28,279265
0,125	0	4,077794902	0,625	0	4,8303102	1,125	0	6,5689854	1,625	0	9,3010573
0,125	0,125	4,083392226	0,625	0,125	4,8345214	1,125	0,125	6,5718746	1,625	0,125	9,3033569
0,125	0,25	4,124655654	0,625	0,25	4,8745571	1,125	0,25	6,6110575	1,625	0,25	9,3426272
0,125	0,375	4,236298249	0,625	0,375	4,9852581	1,125	0,375	6,7213593	1,625	0,375	9,4531903
0,125	0,5	4,453128489	0,625	0,5	5,2014589	1,125	0,5	6,9375673	1,625	0,5	9,6697278
0,125	0,625	4,810075442	0,625	0,625	5,5580652	1,125	0,625	7,2945187	1,625	0,625	10,027092
0,125	0,75	5,342170876	0,625	0,75	6,0900442	1,125	0,75	7,8271524	1,625	0,75	10,560269
0,125	0,875	6,084521561	0,625	0,875	6,8324349	1,125	0,875	8,5704381	1,625	0,875	11,304308
0,125	1	7,07227101	0,625	1	7,8203369	1,125	1	9,5594291	1,625	1	12,294292
0,125	1,125	8,340593862	0,625	1,125	9,0888839	1,125	1,125	10,829217	1,625	1,125	13,56532
0,125	1,25	9,924686593	0,625	1,25	10,673236	1,125	1,25	12,414912	1,625	1,25	15,152501
0,125	1,375	11,85974207	0,625	1,375	12,608557	1,125	1,375	14,351653	1,625	1,375	17,090967
0,125	1,5	14,18095412	0,625	1,5	14,930032	1,125	1,5	16,674587	1,625	1,5	19,415833
0,125	1,625	16,92352093	0,625	1,625	17,672847	1,125	1,625	19,41885	1,625	1,625	22,162225
0,125	1,75	20,12263777	0,625	1,75	20,872196	1,125	1,75	22,619604	1,625	1,75	25,365184
0,125	1,875	23,81349	0,625	1,875	24,563272	1,125	1,875	26,312019	1,625	1,875	29,059835
0,25	0	4,174172984	0,75	0	5,1720918	1,25	0	7,1587763	1,75	0	10,139262
0,25	0,125	4,179452506	0,75	0,125	5,1759439	1,25	0,125	7,1614068	1,75	0,125	10,141851
0,25	0,25	4,220372309	0,75	0,25	5,2157245	1,25	0,25	7,2004768	1,75	0,25	10,181484
0,25	0,375	4,331711049	0,75	0,375	5,3262811	1,25	0,375	7,310766	1,75	0,375	10,292263
0,25	0,5	4,548293067	0,75	0,5	5,5424416	1,25	0,5	7,527037	1,75	0,5	10,508905
0,25	0,625	4,905042335	0,75	0,625	5,8991008	1,25	0,625	7,8841009	1,75	0,625	10,866336
0,25	0,75	5,436976133	0,75	0,75	6,431213	1,25	0,75	8,4169006	1,75	0,75	11,399579
0,25	0,875	6,179189313	0,75	0,875	7,1738005	1,25	0,875	9,1604138	1,75	0,875	12,143708
0,25	1	7,166815957	0,75	1	8,1619466	1,25	1	10,14969	1,75	1	13,133831
0,25	1,125	8,435021644	0,75	1,125	9,4307608	1,25	1,125	11,41982	1,75	1,125	14,405063
0,25	1,25	10,01900509	0,75	1,25	11,015385	1,25	1,25	13,00591	1,75	1,25	15,992541
0,25	1,375	11,95396265	0,75	1,375	12,950982	1,25	1,375	14,943086	1,75	1,375	17,931407
0,25	1,5	14,27508101	0,75	1,5	15,272739	1,25	1,5	17,266479	1,75	1,5	20,256749
0,25	1,625	17,0175547	0,75	1,625	18,015834	1,25	1,625	20,011216	1,75	1,625	23,003709
0,25	1,75	20,21657357	0,75	1,75	21,215457	1,25	1,75	23,212448	1,75	1,75	26,207364
0,25	1,875	23,90732651	0,75	1,875	24,906786	1,25	1,875	26,905329	1,75	1,875	29,902791
0,375	0	4,331570355	0,875	0	5,5757217	1,375	0	7,8107216	1,875	0	11,039132
0,375	0,125	4,336507689	0,875	0,125	5,5792307	1,375	0,125	7,8131474	1,875	0,125	11,042637
0,375	0,25	4,3771121	0,875	0,25	5,6187792	1,375	0,25	7,8521775	1,875	0,25	11,082858
0,375	0,375	4,488197757	0,875	0,375	5,7292175	1,375	0,375	7,9625004	1,875	0,375	11,193907
0,375	0,5	4,704599298	0,875	0,5	5,9453657	1,375	0,5	8,1788491	1,875	0,5	11,410648
0,375	0,625	5,06123568	0,875	0,625	6,3021052	1,375	0,625	8,5360227	1,875	0,625	11,768109
0,375	0,75	5,593106987	0,875	0,75	6,8343777	1,375	0,75	9,0689697	1,875	0,75	12,30136
0,375	0,875	6,335292842	0,875	0,875	7,5771858	1,375	0,875	9,8126882	1,875	0,875	13,045515
0,375	1	7,322919826	0,875	1	8,5656042	1,375	1	10,802235	1,875	1	14,035724
0,375	1,125	8,591141167	0,875	1,125	9,834729	1,375	1,125	12,072698	1,875	1,125	15,307122
0,375	1,25	10,17514142	0,875	1,25	11,419682	1,375	1,25	13,659178	1,875	1,25	16,894837
0,375	1,375	12,11011184	0,875	1,375	13,355619	1,375	1,375	15,596795	1,875	1,375	18,834031
0,375	1,5	14,43124363	0,875	1,5	15,67772	1,375	1,5	17,920665	1,875	1,5	21,15982
0,375	1,625	17,17373441	0,875	1,625	18,421156	1,375	1,625	20,665916	1,875	1,625	23,90734
0,375	1,75	20,37277167	0,875	1,75	21,621108	1,375	1,75	23,867684	1,875	1,75	27,111724
0,375	1,875	24,06354739	0,875	1,875	25,312743	1,375	1,875	27,561086	1,875	1,875	30,808018

 

Полученный график функции  представлен ниже:

 

Проанализируем сходимость решения, рассмотрев конкретную точку  – x=1, y=1.

Число разбиений - Nx	Значение функции – U(x,y)	Разница между значениями
8	9,125380	-
16	9,031412	0,093968
32	9,007793	0,023619
64	9,001894	0,005899

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

 

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

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

k – количество итераций	Nx	Ny	N- порядок матрицы	Погрешность
6	4	4	16	0,780554
10	8	8	64	0,199844
17	16	16	256	0,050320
31	32	32	1024	0,012524
60	64	64	4096	0,003292

 

Вывод

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

 

 

 

Приложение 1

Интерфейс программы

 

 

Входными данными для  программы являются:

• Границы исследуемых диапазонов;
• Количество разбиений;
• Требуемая точность метода;
• Входные функции V1, V2, V3, V4, U, k1, k2, fun;