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

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

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

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


 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Отчет

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

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

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

 

 

 

 

 

 

 

 

 

Выполнила:

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

Гр. 3084/3

Проверил:

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

 

 

Содержание

 

 

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

 

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

 

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

 

          

 

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

 

 




 




 

 


 

 

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

 

 

 

 

 

 

 

 

 

 

 

 

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

 

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

 

a=xL<x1< ...<xR=b

c=yL<y1< ... <yR=d

 

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

                                                          (2)

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

   

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

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

 

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

 

 (1)

 

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

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

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

(2)

 

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

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

  (3)   

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

 

Для x=x0 y=y0:

 

i=0, j=0

 

 (5)

 

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

 

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

 

Vi,j = Wk

Vi-1,j = Wk-1

Vi+1,j = Wk+1

Vi,j-1  = Wk-L

Vi,j+1 = Wk+L,  где L = Nx

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

 

k =  j*Nx + i,

 

где i = 0..Nx-1;   j = 0..Ny-1.

 

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

N=Nx*Ny.

 

 

 

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

 

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

 

A*X=F

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

 

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

 

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

 

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

 




 

      А =


 

 

 

 

 

 

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

 

 

 

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

 

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

N=Nx*Ny

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

 

i=1..Nx-1,  j=1..Ny-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;

g1(y) = 3;

g2(y) = 3;

g3(x) = 3;

g4(x) = 3;

k1(x,y) = 2; k2(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;

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

g2(y) = (2*xRight+3*y+5);

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

g4(x) = 2*x+3*yRight+5;

k1(x,y) = 2*x + y + 1;

k2(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;

k1(x,y) = 1 + 2*x*x + y;

k2(x,y) = 1 + x + 2*y*y;

; ;

 

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

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

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

 g2(y) = 2*xRight*xRight+3*y*y*y+4;

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

g4(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;

Также дополнительно  можно установить режим анализа  погрешности, позволяющий проследить изменение погрешности. Результаты по этому анализу заносятся в таблицу на вкладке “Анализ погрешности”:

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

На этой вкладке также  выводится дополнительная информация, такая как

  • максимальная погрешность
  • количество итераций метода
  • число разбиений Nx, Ny и N
  • потраченное время

Текст программы

Определение типов данных и методов

 

using System;

using System.IO;

namespace MMKursProject

{

public class Matrix

{

public int N;

public int Nx;

public bool fiveDiag;

public bool econom;

public bool transpon;

public double[] A, B, C;

public double[][] M;

 

public Matrix(int N, int Nx, bool econom, bool fiveDiag)

{

this.Nx = Nx;

this.N=N;

this.fiveDiag=fiveDiag;

this.econom=econom;

transpon=false;

if(econom)

{ A = new double[N];

B = new double[N - 1];

C = new double[N - Nx];

}

else

{ M=new double[N][];

for(int i=0;i<N;i++)

M[i]=new double[N];

}

}

public Matrix(int N, int Nx, string str)

{

this.Nx = Nx;

this.N=N;

transpon=false;

if(str=="econom") this.econom=true;

else if(str=="full") this.econom=false;

else throw new ArgumentOutOfRangeException();

if(econom)

{

A = new double[N];

B = new double[N - 1];

C = new double[N - Nx];

}

else

{

M=new double[N][];

for(int i=0;i<N;i++)

M[i]=new double[N];

}

}

public Matrix(int N, int Nx, double[] A, double[] B, double[] C)

{

if ((A.Length!=N)||(B.Length!=N-1)||(C.Length!=N-Nx))

throw new ArgumentOutOfRangeException();

this.Nx = Nx;

transpon=false;

this.N=N;

this.econom=true;

this.A = new double[N];

this.B = new double[N - 1];

this.C = new double[N - Nx];

for (int i=0;i<N;i++)  this.A[i]=A[i];

for (int i=0;i<N-1;i++) this.B[i]=B[i];

for (int i=0;i<N-Nx;i++) this.C[i]=C[i];

}

public Matrix(int N, int Nx, Vector A, Vector B, Vector C)

{

if ((A.GetLenght!=N)||(B.GetLenght!=N-1)||(C.GetLenght!=N-Nx))

throw new ArgumentOutOfRangeException();

this.Nx = Nx;

this.N=N;

transpon=false;

this.econom=true;

this.A = new double[N];

this.B = new double[N - 1];

this.C = new double[N - Nx];

for (int i=0;i<N;i++)  this.A[i]=A[i];

for (int i=0;i<N-1;i++) this.B[i]=B[i];

for (int i=0;i<N-Nx;i++) this.C[i]=C[i];

}

public int GetLength

{

get { return A.Length; }

}

public Matrix(StreamReader sr): this(Int16.Parse(sr.ReadLine()),Int16.Parse(sr.ReadLine()),sr.ReadLine())

{ string str;

transpon=false;

int i=0,j=0,k=0,newk=0;

if(!econom)

{ while (((str=sr.ReadLine())!=null)&&(i<N))

{ j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf('\t',k+1);

this[i,j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while(k!=-1);

i++;

}

}

else

{ fiveDiag = true;

if ((str=sr.ReadLine())!=null)

{

j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

A[j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<N));

}

else throw new ArrayTypeMismatchException();

if ((str=sr.ReadLine())!=null)

{

j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

B[j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<N-1));

}

else throw new ArrayTypeMismatchException();

if ((str=sr.ReadLine())!=null)

{

j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

C[j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<N-Nx));

}

else throw new ArrayTypeMismatchException();

}

}

 

public double this[int i, int j]

{

get {return GetElement(i, j);}

set {SetElement(i, j, value);}

}

 

public double GetElement(int i, int j)

{ if(transpon)

{ int temp=i;

i=j;

j=temp;

}

if ((i < 0) || (j < 0) || (i >= N) || (j >= N))// Выход за границы

throw new ArgumentOutOfRangeException();

if(econom)

{ if (fiveDiag) // если 5-диагональная матрица

{ if (i==j)// Диагональ

return A[i];

if(i>j)// Нижний треугольник

{

if(i==j+1) return B[j];

if(i==j+Nx) return C[j];

}

else// Верхний треугольник

{

if(i==j-1) return B[i];

if(i==j-Nx) return C[i];

}

}

else //нижняя треугольная трехдиагональная матрица

{

if (i==j)// Диагональ

return A[i];

// Нижний треугольник

if (i > j)

{

if(i==j+1) return B[j];

if(i==j+Nx) return C[j];

}

 

}

return 0;

}

else return M[i][j];

}

public void SetElement(int i, int j, double value)

{ if(transpon)

{ int temp=i;

i=j;

j=temp;

}

if ((i < 0) || (j < 0) || (i >= N) || (j >= N))// Выход за границы

throw new ArgumentOutOfRangeException();

if(econom)

{ if (fiveDiag)

{ if (i==j)// Диагональ

A[i]=value;

if(i>j)// Нижний треугольник

{

if(i==j+1) B[j]=value;

if(i==j+Nx) C[j]=value;

}

else// Верхний треугольник

{

if(i==j-1) B[i]=value;

if(i==j-Nx) C[i]=value;

}

}

else

{ if (i==j)// Диагональ

A[i]=value;

// Нижний треугольник

 

if((i==j+1)&&(i>j)) B[j]=value;

if((i==j+Nx)&&(i>j)) C[j]=value;

if((i==j-1)&&(i<j)) B[i]=value;

if((i==j-Nx)&&(i<j)) C[i]=value;

 

}

}

else M[i][j]=value;

}

public static Matrix operator *(Matrix m1, Matrix m2)

// умножение матрицы  на матрицу

{ Matrix res = null;

int len=0;

if ((len=m1.N) != m2.N) throw new ArgumentOutOfRangeException();

else

{ res = new Matrix(len,m1.Nx,false,true);

for (int n = 0; n < len; n++)

{ for (int i=0; i<len; i++)

{ res[n,i]=0;

for (int j = 0; j < len; j++)

{res[n,i] += m1[n, j]*m2[j,i];}

}

}

}

return res;

}

public override string ToString()

{ string str="";

for(int i=0; i<N;i++)

{for(int j=0; j<N; j++) str+=this[i,j].ToString("#,##0.00")+

((j<N-1)?"\t":"");//"#,###0.000"

str+='\n';

}

return str;

}

}

public class Vector

{ public double[] V;

public Vector(int size)

{

V=new double[size];

}

public Vector(int size, double[] V)

{

this.V=new double[size];

for(int i=0;i<size;i++)

this.V[i]=V[i];

}

public int GetLenght

{ get { return V.Length; }

}

public Vector(StreamReader sr):this(Int16.Parse(sr.ReadLine()))

{ string str;

int j=0,k=0,newk=0;

if ((str=sr.ReadLine())!=null)

{ j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

V[j]=double.Parse(newk>0?str.Substring(k,newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<V.Length));

}

else throw new ArrayTypeMismatchException();

}

 

public double this[int i]

{ get

{ if(i>=V.Length) throw new ArgumentOutOfRangeException();

return V[i];

}

set

{ if(i>=V.Length) throw new ArgumentOutOfRangeException();

V[i]=value;

}

}

public static double operator *(Vector v1, Vector v2)

//умножение вектора  на вектор

{ double res = 0;

int len=0;

if ((len=v1.V.Length) != v2.V.Length) throw new ArgumentOutOfRangeException();

else

for (int i = 0; i < len; i++)

{ res += v1[i]*v2[i];

}

return res;

}

public static Vector operator +(Vector v1, Vector v2)

// сумма двух векторов

{ Vector res = null;

int len=0;

if ((len=v1.V.Length) != v2.V.Length) throw new ArgumentOutOfRangeException();

else

{ res = new Vector(len);

for (int i = 0; i < len; i++)

{res[i] = v1[i] + v2[i];

}

}

return res;

}

public static Vector operator -(Vector v1, Vector v2)

// разность двух  векторов

{

Vector res = null;

int len=0;

if ((len=v1.V.Length) != v2.V.Length) throw new ArgumentOutOfRangeException();

else

{

res = new Vector(len);

for (int i = 0; i < len; i++)

{

res[i] = v1[i] - v2[i];

}

}

return res;

}

public static Vector operator *(Matrix m, Vector v)

// умножение матрицы  на вектор

{

Vector res = null;

int len=0;

if ((len=m.N) != v.V.Length) throw new ArgumentOutOfRangeException();

else

{

res = new Vector(len);

for (int i = 0; i < len; i++)

{

res[i] = 0;

for (int j = 0; j < len; j++)

{res[i] += m[i, j]*v[j];

}

}

}

return res;

}

public static Vector operator *(Vector v, Matrix m)

// умножение вектора  на матрицу

{ Vector res = null;

int len=0;

if ((len=m.N) != v.V.Length) throw new ArgumentOutOfRangeException();

else

{

res = new Vector(len);

for (int j = 0; j < len; j++)

{ res[j] = 0;

for (int i = 0; i < len; i++)

{res[j] += m[i, j]*v[i];

}

}

}

return res;

}

public static Vector operator *(double d, Vector v)

//умножение скаляра  на вектор

{ int len=v.V.Length;

Vector res = new Vector(len);

for (int i = 0; i < len; i++)

{res[i] = v[i]*d;

}

return res;

}

public static Vector operator *(Vector v, double d)

//умножение вектора  на скаляр

{ int len=v.V.Length;

Vector res = new Vector(len);

for (int i = 0; i < len; i++)

{res[i] = v[i]*d;

}

return res;

}

public static double Abs(Vector v)

{return (Math.Sqrt((v*v)));

}

 

public override string ToString()

{ string str="";

int len=V.Length;

for(int i=0;i<len;i++) str+=V[i].ToString("#,#####0.00000")+" ";

return str;

}

}

}

 

Метод сопряженных градиентов

 

using System;

using System.Windows.Forms;

namespace MMKursProject

{

public class Solving

{ public static Vector X;

public static int Nx = 4;

public static int Ny = 4;

public static double Eps=0.00000001;

public static double xLeft=0, xRight=1, yLeft=0, yRight=1;

public static double hx, hy;

public static int N,i,j;

public static double kappa1 = 1;

public static double kappa2 = 1;

public static double kappa3 = 1;

public static double kappa4 = 1;

public static double[] xi, yj;

public static double[] Ck, Bk, Ak, Fk,Uk,XX;

public static int k = 0, iter = 0;

public static double maxerror = 0;

public static DateTime start;

public static TimeSpan ts;

public static int State=0;

 

public static void Solve()

{

try

{

start = DateTime.Now;

N = Nx*Ny;

hx = (xRight - xLeft)/Nx;

hy = (yRight - yLeft)/Ny;

xi = new double[Nx + 1];

yj = new double[Ny + 1];

// сетка по х

for (i = 0; i <= Nx; i++)

{ xi[i] = xLeft + i*hx;  }

// сетка по у

for (j = 0; j <= Ny; j++)

{ yj[j] = yLeft + j*hy;  }

Ck = new double[N];

Bk = new double[N];

Ak = new double[N];

Fk = new double[N];

Uk = new double[N];

XX = new double[N];

for (j=0; j<=Ny-1; j++)

for (i=0; i<=Nx-1; i++)

{

k=j*Nx+i;

Uk[k]=u(xi[i],yj[j]);

}

 

Ak[0]=(hy/2)*(k1(xi[0],xi[1],yj[0])/hx + kappa1) + (hx/2)*(k2(xi[0],yj[0],yj[1])/hy + kappa3);

Fk[0]=(hx*hy/4)*fun(xi[0],yj[0]) + (hy/2)*V1(yj[0]) + (hx/2)*V3(xi[0]);

for ( i = 1; i<=(Nx-1); i++)

for (j = 1; j<= (Ny-1); j++)

{

k = j*Nx + i;

Ak[k]= (hy/hx)*(k1(xi[i],xi[i+1],yj[j]) + k1(xi[i],xi[i-1],yj[j])) + (hx/hy)*(k2(xi[i],yj[j],yj[j+1]) + k2(xi[i],yj[j],yj[j-1]));

Bk[k]= -(hy/hx)*k1(xi[i],xi[i-1],yj[j]) ;

Ck[k]= -(hx/hy)*k2(xi[i],yj[j],yj[j-1]);

Fk[k]= hx*hy*fun(xi[i],yj[j]);

if (i==(Nx-1)) Fk[k]= Fk[k] + (hy/hx)*V2(yj[j])*k1(xi[i],xi[i+1],yj[j]);

if (j==(Ny-1)) Fk[k]= Fk[k] + (hx/hy)*V4(xi[i])*k2(xi[i],yj[j],yj[j+1]); 

}

for (j = 1; j<=(Ny-1); j++)

{

i = 0;

k=j*Nx + i;

Ak[k]=hy*(k1(xi[i],xi[i+1],yj[j])/hx + kappa1) + (hx/(2*hy))*(k2(xi[i],yj[j],yj[j+1]) + k2(xi[i],yj[j],yj[j-1]));

Ck[k]= -(hx/(2*hy))*k2(xi[i],yj[j],yj[j-1]);

Fk[k]=(hx/2)*hy*fun(xi[i],yj[j]) + hy*V1(yj[j]);

if (j==(Ny-1)) Fk[k]= Fk[k] + (hx/(2*hy))*V4(xi[i])*k2(xi[i],yj[j],yj[j+1]);

}

for (i = 1; i<=(Nx-1); i++)

{

j = 0;

k=j*Nx + i;

Ak[k]=hx*(k2(xi[i],yj[j],yj[j+1])/hy + kappa3) + (hy/(2*hx))*(k1(xi[i],xi[i+1],yj[j]) + k1(xi[i],xi[i-1],yj[j]));

Bk[k]= -(hy/(2*hx))*k1(xi[i],xi[i-1],yj[j]);

Fk[k]=hx*(hy/2)*fun(xi[i],yj[j]) + hx*V3(xi[i]);

if (i==(Nx-1)) Fk[k]= Fk[k] + (hy/(hx*2))*V2(yj[j])*k1(xi[i],xi[i+1],yj[j]);

}

 

//------------------------------------------------------------------------------------

// формируем вектора  для матрицы

Vector A = new Vector(N);

Vector B = new Vector(N - 1);

Vector C = new Vector(N - Nx);

Vector F = new Vector(N);

Vector X0 = new Vector(N);

Vector UU = new Vector(N);

for (k = 0; k < N; k++) A[k] = Ak[k];

for (k = 0; k < N - 1; k++) B[k] = Bk[k + 1];

for (k = 0; k < N - Nx; k++) C[k] = Ck[k + Nx];

for (k = 0; k < N; k++) F[k] = Fk[k];

for (k = 0; k < N; k++) UU[k] = Uk[k]; 

// формируем матрицу

Matrix MM = new Matrix(N, Nx, A, B, C);

MM.fiveDiag = true;

//X = Metod(MM, F, X0); // метод сопряженных градиентов без предобусловливания

X = MetodHoleckogo(MM, F, X0); // метод сопряженных градиентов без предобусловливания

maxerror=Math.Abs(UU[0]-X[0]);

for (i=1;i<N;i++)

{

if(maxerror<Math.Abs(UU[i]-X[i]))

{ maxerror=Math.Abs(UU[i]-X[i]); }

}

DateTime end = DateTime.Now;

ts=end-start;

State=2;

 

}

catch (Exception exc)

{ MessageBox.Show(exc.Message+"***");  }

}

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

public static Vector Metod(Matrix A, Vector F, Vector x0)

{ k = 0;

bool Stop = false;

Vector r0 = new Vector(F.GetLenght);

Vector r1 = new Vector(F.GetLenght);

Vector S1 = new Vector(F.GetLenght);

Vector x = new Vector(F.GetLenght);

Vector x00 = new Vector(F.GetLenght);

Double alfa, betta;

x00 = x0;

k = 0;

// невязка

r0 = (F - A*x0);

//  вспомогательный  вектор в методе

S1 = r0;

while (!Stop)

{

k++;

// подсчет альфы

alfa = (r0*r0)/((A*S1)*S1);

x = x0 + alfa*S1;

r1 = r0 - alfa*(A*S1);

// проверка на  сходимость - завершение метода

if ((Vector.Abs(r1)/Vector.Abs(F)) < 0.000001)

{ if ((Vector.Abs(F - A*x00)/Vector.Abs(F)) > 0.000001)

{ Stop = true;  }

else return x;

}

betta = (r1*r1)/(r0*r0);

S1 = r1 + betta*S1;

r0 = r1;

x0 = x;

}

return x;

}

// метод сопряженных градиентов с предобуславливанием

public static Vector MetodHoleckogo(Matrix A, Vector F, Vector x0)

{

iter = 0;

bool Stop = false;

Vector r0 = new Vector(F.GetLenght);

Vector r1 = new Vector(F.GetLenght);

Vector S1 = new Vector(F.GetLenght);

Vector x = new Vector(F.GetLenght);

Vector x00 = new Vector(F.GetLenght);

Vector w0 = new Vector(F.GetLenght);

Vector w1 = new Vector(F.GetLenght);

Matrix L;

double alfa, betta;

x00 = x0;

iter = 0;

// невязка

r0 = (F - A*x0);

L = CreateMatrixL(A);

w0 = FindW(L, r0);

//  вспомогательный  вектор в методе

S1 = w0;

while (!Stop)

{

iter++;

// подсчет альфы

alfa = (w0*r0)/((A*S1)*S1);

x = x0 + alfa*S1;

r1 = r0 - alfa*(A*S1);

// проверка на  сходимость - завершение метода

if ((Vector.Abs(r1)/Vector.Abs(F)) < 0.000001)

{

if ((Vector.Abs(F - A*x00)/Vector.Abs(F)) > 0.000001)

{Stop = true;}

else return x;

}

w1 = FindW(L, r1);

betta = (w1*r1)/(w0*r0);

S1 = w1 + betta*S1;

r0 = r1;

w0 = w1;

x0 = x;

}

return x;

}

public static Matrix CreateMatrixL(Matrix M)

{ //Получаем вектора a,b,c из основной матрицы для формирования векторов в матрицу L

Vector A = new Vector(M.GetLength);

Vector B = new Vector(M.GetLength);

Vector C = new Vector(M.GetLength);

double[] Al = new double[M.GetLength];

double[] Bl = new double[M.GetLength];

double[] Cl = new double[M.GetLength];

for (int i = 0; i < M.N; i++)

{ for (int j = 0; j < M.N; j++)

{ if (i == j) // Диагональ

A[i] = M[i, j];

if (i > j) // Нижний треугольник

{

if (i == j + 1) B[j] = M[i, j];

if (i == j + M.Nx) C[j] = M[i, j];

}

}

}

Al[0]=Math.Sqrt(A[0]);

Bl[0]=B[0]/Al[0];

Cl[0] = C[0]/Al[0];

for (int i = 1; i < M.Nx; i++)

{ Al[i] = Math.Sqrt(A[i]-Bl[i-1]*Bl[i-1]);

Bl[i] = B[i]/Al[i];

Cl[i] = C[i]/Al[i];

}

for (int i = M.Nx; i < M.N; i++)

{ Al[i] = Math.Sqrt(A[i]-Bl[i-1]*Bl[i-1]-Cl[i-M.Nx]*Cl[i-M.Nx]);

Bl[i] = B[i]/Al[i];

Cl[i] = C[i]/Al[i];

}

//формируем вектора  для построения матрицы

A = new Vector(M.N);

B = new Vector(M.N - 1);

C = new Vector(M.N - M.Nx);

for (int i = 0; i < M.N; i++) A[i] = Al[i];

for (int i = 0; i < M.N - 1; i++) B[i] = Bl[i];

for (int i = 0; i < M.N - M.Nx; i++) C[i] = Cl[i];

Matrix L = new Matrix(M.N, M.Nx, A, B, C);

return L;

}

// метод  Гаусса - решает систему вниз

public static Vector GaussDown(Matrix M, Vector F)

{ Vector x = new Vector(M.N);

x[0] = F[0]/M.GetElement(0, 0);

for (int i = 1; i < M.Nx; i++) x[i] = (F[i] - M.GetElement(i, i - 1)*x[i - 1])/M.GetElement(i, i);

for (int i = M.Nx; i < M.N; i++)

x[i] = (F[i]-M.GetElement(i, i-M.Nx)*x[i-M.Nx]-M.GetElement(i,i-1)*x[i-1])/M.GetElement(i, i);

return x;

}

// метод  Гаусса - решает систему вверх

public static Vector GaussUp(Matrix M, Vector F)

{ M.transpon = true;

Vector x = new Vector(M.N);

x[M.N - 1] = F[M.N - 1]/M.GetElement(M.N - 1, M.N - 1);

for (int i = M.N - 2; i >= (M.N - M.Nx); i--)

x[i] = (F[i] - M.GetElement(i, i + 1)*x[i + 1])/M.GetElement(i, i);

for (int i = M.N - M.Nx - 1; i >= 0; i--)

x[i] = (F[i]-M.GetElement(i,i+M.Nx)*x[i+M.Nx]-M.GetElement(i,i+1)*x[i+1])/M.GetElement(i,i);

M.transpon = false;

return x;

}

// нахождение дополнительного вектора w

public static Vector FindW(Matrix L, Vector r)

{ Vector temp = new Vector(L.GetLength, GaussDown(L, r).V);

return GaussUp(L, temp);

}

// функция к1

public static double k1(double xi, double xi_1, double yi)

{ double x = (xi + xi_1)/2;

double res;

//res = 2; //test 1

//res = 2*x + yi + 1; //test 2

res = 1 + 2*x*x + yi; //test3

return res;

}

// функция к2

public static double k2(double xi, double yi, double yi_1)

{ double res;

double y = (yi + yi_1)/2;

//res = 5; //test 1

//res = xi + 2*y + 1; //test2;

res = 1 + xi + 2*y*y; //test3;

return res;

}

//решение уравнения - u

public static double u(double x, double y)

{

double res;

//res = 3; // test1

//res = 2*x+3*y+5; //test2

res = 2*x*x + 3*y*y*y + 4; //test3;

return res;

}

 

public static double V1(double y)

{

double res;

//res = 3; // test1

//res =kappa1*(2*xLeft+3*y+5) - (2*xLeft+y+1)*2; //test2

res = kappa1*(2*xLeft*xLeft+3*y*y*y+4)-(1+2*xLeft*xLeft+y)*(4*xLeft);//test3

return res;

}

 

public static double V2(double y)

{

double res;

Математические модели технических объектов 📙 Курсовая → 🆔 121098

Математические модели технических объектов 📙 Курсовая → 🆔 121098