Рома_Шелев_II

Содержание

5.5вспомогательная программа GetA. c . . . . . . . . . . . . . . . . . . . . . . . 12

5.6Функции на Matlab для построения графиков . . . . . . . . . . . . . . . . . . 12

2

1Постановка учебно-практической задачи

Решить для n=10000 методом прогонки и методом минимальной невязки систему уравнений с относительной точностью 0.001.

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

2Описание алгоритма

2.1Квадратурная формула трапеций (n=1)

2.2Рекуррентная квадратурная формула трапеций

где h = 2−N и точки xk = α + kh делят интервал [α, β] на 2N = 2M точек.

2.3Метод прогонки

Важным классом систем являются системы вида:

d1x1 + c1x2 = b1

aixi−1 + dixi + cixi+1 = bi, i = 2, n 1 anxn−1 + dnxn = bn

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

3

x3. Следовательно третье уравнение дает соотношение между x3 и x4 и т.д. Запишем связь

Вычисляем коэффициенты Li, Mi(i = 2, n), процесс нахождения которых называется прямым ходом метода прогонки.

1 + sin(i)2 Li+1 = Li (i2 + 5)

bi Mi

Mi+1 = Li (i2 + 5), i = 2, n 1

∫ inf

bi = (cos(2t)exp( i2 t2 t2))dt

0

2.4Метод минимальной невязки

Вданном методе:

Hk = τkE, Tk = E τkA,

В методе минимальных невязок итерационный параметр τk выбирается из условия минимума нормы невязки на (k + 1)-ом шаге и имеет вид

4

2.5Метод простой итерации

Пусть система линейных алгебраических уравнений (СЛАУ) задана в виде:

стемы, разрешенной относительно i-й перемнной. Рабочие формулы метода простой итерации имеют вид:

∫ inf

bi = (cos(2t)exp( i2 t2 t2))dt

0

3Исследование применимости алгоритма

3.1Квадратурная формула трапеций

Исследуем применимость квадратурной формулы трапеций к интегралу

∫ inf

(cos(2t)exp( i2 t2 t2))dt

0

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

fn(x) ! f(x), при n ! 1,

т. е. при увеличении числа узлов сетки ωn значение интеграла будет стремиться к реальному значению. Соответственно, для каждого значения i можно оценить значение остаточного члена.

Точность вычисления контролировалась по правилу Рунге (т.е. сравнивались значения, полученные при сетке с шагом h и h/2).

5

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

с относительной точностью 0.01. Для обеспечения необходимой точности необходимо выполнение неравенства:

jjbk bk+1jj < ϵ,

bk+1

где ϵ - относительная точность вычисления интеграла (в нашем случае 0.01).

3.2Метод простой итерации

Рассмотрим применимость метода наискорейшего к исследуемой системе уравнений. С помощью итерационного метода ищется приближенное решение системы. Говорят,

что итерация xk является с точностью ϵ приближенным решением системы, если

jjxk x0jj ϵ

где x0 — точное решение системы. При этом норма jjxjj вектора x = fx1, ..., xng 2 Rn равна jjxjj∞ = max1≤i≤n jxij, ϵ - абсолютная погрешность.

Итерационный метод называется стационарным, если матрица Hk не зависит от номера шага k. В противном случае метод называется нестационарным.Для того, чтобы стационарный итерационный процесс сходился, достаточно, чтобы для какой-либо одной нормы матрицы T выполнялось неравенство

jjT jj < 1.

Исследуем сходимость системы линейных уравнений. Можно оценить норму матрицы T по формуле

При i = 1, n и j = 1, n . В каждой строке матрицы Т находится два элемента с коэф-

max i + 4

1≤i≤n i2 + 5

построив график функции y = xx2+4+5 мы видим, что максимальное значение он принимает на x = 0.5 равное 0,86, а на бесконечности стремиться к 0.

отсюда следует, что jjT jj < 1, что говорит о сходимости алгоритма.

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

jjxk xk−1jj < ϵ, jjxkjj

где ϵ - относительная точность решения системы линейных уравнений.

4Результаты расчетов

4.1Графики

6

4

3.5

3

2.5

2

1.5

1

0.5

0

Рис. 2: Подсчет интегралов при n = 10000

7

4

3.5

3

2.5

2

1.5

1

0.5

0

Рис. 4: Решение линейной системы при n = 10000

8

5Приложение

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

5.1Вычисление интеграла по формуле трапеции

9

/ ( ( ( double ) ( I I ) + 5 . ) ) ;

} ;

5.2Метод минимальной невязки

include <stdio.h> include <math.h> include "Trapezy.c"include "GetA.c" define N 20 define Accuracy 0.01

double Delta ; double Tau = 1. ; // Tau parametr of system double X[ N ] ; // Solution of system double R[ N ] ; // R array double AR[ N ] ;// A*R array double B[ N ] ; // Integral array

void CalcR() int I , J ;

R[ 0 ] = GetA( 1 , 1 ) * X[ 0 ] + GetA( 1 , 2 ) * X[ 1 ] - B[ 0 ] ; R[ N - 1 ] = GetA( N , N ) * X[ N - 1 ] - B[ N - 1 ] ; for ( I = 1 ; I < N - 1 ; I++ ) R[ I ] = 0. ; for ( J = I - 1 ; J < I + 2 ; J++ ) R[ I ] += GetA( I + 1 , J + 1 ) * X[ J ] ; R[ I ] -= B[ I ] ; ; ;

void CalcAR() int I , J ;

AR[ 0 ] = GetA( 1 , 1 ) * R[ 0 ] + GetA( 1 , 2 ) * R[ 1 ] ; AR[ N - 1 ] = GetA( N , N ) * R[ N - 1 ] ; for ( I = 1 ; I < N - 1 ; I++ ) AR[ I ] = 0. ; for ( J = I - 1 ; J < I + 2 ; J++ ) AR[ I ] += GetA( I + 1 , J + 1 ) * R[ J ] ; ;

;

void CalcTau() int I ; double Ch = 0. , Zn = 0. ;

for ( I = 0 ; I < N ; I++ ) Ch += AR[ I ] * R[ I ] ; Zn += AR[ I ] * AR[ I ] ; ; Tau = Ch / Zn ; ;

void CalcX() double NewX[ N ] ; int I , J ; NewX[ 0 ] = ( 1. - Tau * GetA( 1 , 1 ) ) * X[ 0 ] - Tau * GetA( 1 , 2 ) * X[ 1 ] + Tau * B[ 0 ] ; NewX[ N - 1 ] = ( 1. - Tau * GetA( N , N ) ) * X[ N - 1 ] + Tau * B[ N - 1 ] ; for ( I = 1 ; I < N - 1 ; I++ ) NewX[ I ] = -Tau * GetA( I + 1 , I ) * X[ I - 1 ] + ( 1. - Tau * GetA( I + 1 , I + 1 ) ) * X[ I ] - Tau * GetA( I + 1 , I + 2 )

/* Functions to out Solution */ /* ======================= */ void PrintX() int I ; for ( I = 0 ; I < N ; I++ ) printf( "

int main() int I = 0 ;

CalcB( B ) ; for ( I = 0 ; I < N ; I++ ) X[ I ] = 1. ;

do CalcR() ; CalcAR() ; CalcTau() ; CalcX() ; while ( Delta > 0.001 ) ; PrintX() ; return 0 ; ;

5.3Метод прогонки

double Delta ;

10