metodichka_2_semestr

AX B

или

a11x1 a12 x2 ...	a1n xn b1;
a21x1 a22 x2 ...	a2n xn b2 ;
................................................
an1x1 an2 x2 ...	ann xn bn ;	(1)
Разрешим первое уравнение системы (1) относительно		x1 , второе –
относительно x2 и т. д.

x1 1 / a11(b1 0x1 a12 x2 a13 x3 ... a1n xn );

x2 1 / a22 (b2 a21x1 0x 2 a23x3 ... a2n xn );

............................................................................

xn 1 / ann (bn an1x1 an2 x2 ... an,n 1xn 1 0xn ).

или в матричной форме

				X CX D									(2)
		0	a12 / a11	a13 / a11 .....		a1 n / a11					b1 / a11
											b2 / a22
	C	a21 / a22	0	a23 / a22 .....		a2 n / a22			D
		.....	.....	.....	.....	.....					.......
		an1 / ann	an2 / ann	an3 / ann .....		0		;			bn / ann
									x(0)	, x(0) , x		(0)	,...., x(0).
	Зададимся начальными приближениями корней системы (2)
									1	2		3	n

2

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

приближений в правые части системы (2) X (1) CX (0) D. Найденные первые

приближения можно использовать для нахождения вторых и т. д. Итерационный

процесс продолжается до тех пор , пока все			x(k )
			i не станут достаточно близки к
x(k 1)	,i 1,2,3,.....,n,	т.е.
i

M k max xi(k ) xi(k 1)

Для выполнения условий сходимости итерационного процесса необходимо,

чтобы G 1 , т.е. сумма модулей отношений коэффициентов любой строки к диагональному была меньше 1

n

aij

n

1 или

aii

aij

aii

, i 1,2,....n

(3)

j 1, j i

j 1, j i

Если коэффициенты решаемой системы не удовлетворяют условию сходимости,

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

Алгоритм решения системы линейных алгебраических уравнений методом простых итераций при выполнении лабораторной работы записывается следующим образом:

1. Вводятся исходные данные задачи – порядок системы n , точность

2

нахождения корней EPS , значения элементов массивов B и A ,

удовлетворяющих условию (3), значение максимального количества итераций

MM.

2.Вычисляются элементы массивов C и D по формуле (2), при этом

сохраняются массивы B и A , которые затем используются для проверки решения с помощью стандартной программы simq. Задаётся начальное

приближенное решение X (например, xi di ) .

3.Задаётся начальное значение счётчика итераций M 0 .

4.Вычисляется последующее приближение Y CX D .

5.Изменяется содержимое счётчика итераций на 1: M M 1.

6.Выполняется проверка на сходимость метода M MM ; если условие выполняется, то печатается соответствующее сообщение и вычисления прекращаются.

7.Вычисляется максимальная погрешность текущего приближения по формуле

E max Yi Xi .

8.Задаются значения для выполнения последующего приближения по формуле

Xi Yi .

9.Выполняется проверка на достижение заданной точности определения корней

2

E ; если точность не достигнута, то итерации продолжаются с пункта 4.

10.Исходная система уравнений AX B решается с помощью стандартной подпрограммы SIMQ.

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

Программу на Фортране составить из головной программы и нескольких подпрограмм.

Подпрограммы должны содержать: вычисление элементов массивов C и D и

задание начального приближения X (пункт алгоритма 2);		перемножение
	n
матрицы на вектор и сложение результата с вектором	yi Cij	xj di , i 1,.., n
	j 1

(пункт 4); вычисление максимального значения модуля разности элементов двух одномерных массивов (пункт 6).

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

subroutine izm(a,b,c,d,x,n) real a(n,n),b(n),c(n,n),d(n),x(n) integer n,i,j

do i=1,n do j=1,n

c(i,j)=-a(i,j)/a(i,i) end do

c(i,i)=0.

d(i)=b(i)/a(i,i)

x(i)=d(i) end do

2

return

end subroutine izm

subroutine umn(c,x,d,y,n) real c(n,n),d(n),x(n),y(n) integer n,i,j

do i=1,n y(i)=0.

do j=1,n y(i)=y(i)+c(i,j)*x(j) end do

y(i)=y(i)+d(i) end do

return

end subroutine umn

function emax(y,x,n) real y(n),x(n),e,e1,emax integer n,i e=abs(y(1)-x(1))

do i=2,n e1=abs(y(i)-x(i)) if(e<e1)then e=e1

end if end do emax=e

end function emax

program lab5

integer, parameter::n=4

real a(n,n),b(n),c(n,n),d(n),x(n),y(n) real e,eps,emax

integer m,mm,ier,i,j open(1,file='lab5.txt') read(1,1)eps,((a(i,j),j=1,n),i=1,n),b,mm

2

1 format(4x,f8.4/4(/4f8.4)//4f8.4/3x,i2) call izm(a,b,c,d,x,n)

m=0

e=1.

do while (e>eps) call umn(c,x,d,y,n) m=m+1

if(m>mm) then print 3

3 format('Метод расходится') stop

end if e=emax(y,x,n) do i=1,n x(i)=y(i)

end do end do

call simq(a,b,n,ier) write (1,6)m,x,b print 6,m,x,b

6 format(1x,'Количество итераций=',i2/1x,'Приближенное& значение'/4e15.4/1x,'Решение по методу Гаусса'/4e15.4) end program lab5

Исходные данные задачи помещаются в файл Lab5.txt в соответствии с оператором FORMAT c меткой 1.

Пример записи исходных данных в файл приведен ниже.

ESP=0.0001

матрица А

11.,0.,-1.,1.

1.,7.,2.,-1.

2

y(xi 1)

1.,1.,8.,3.

2.,-3.,1.,10.

Вектор В

31.,3.,-1.,10.

ММ 50

Результаты счета также выводятся в файл Lab5.txt

Лабораторная работа №6

Решение задачи Коши

Задача Коши формулируется следующим образом: на отрезке x0 ; xk

требуется найти решение y y(x) обыкновенного дифференциального

dx f (x, y)	, удовлетворяющее начальному условию y0 (x0 ) y0
уравнения dy

. Численное решение задачи состоит в построении таблицы

	приближенных значений y1, y2 , y3 ,..., y m решения уравнения				y(x) в точках
	x1, x2 ,..., x m .	Обычно	принимают xi x0 i h ,	i 1,2,...,m .	Точки xi
					h	(xk x0 )
	называются	узлами	сетки, а величина	h - шагом:		m .

Существуют два широких класса методов приближенного решения задачи Коши.

1.Одноступенчатые (одношаговые) методы, в которых для

нахождения следующей точки требуется информация лишь об

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

2

y(x)

точности.

2.Многоступенчатые (многошаговые) методы, в которых для отыскания следующей точки решения y(xi 1) используется информация более чем об одной из предыдущих точек и для достижения достаточной точности требуются итерации. Большинство методов этого класса называются методами прогноза и коррекции (методы Милна, Адамса и др.). Далее рассматриваются только одношаговые методы.

Разложение решения y y(x) в окрестности точки x в ряд Тейлора имеет вид

yi 1 yi h yi ' h2 yi ''' h3 yi ''' ..., 2! 3!

где yi 1 y (xi ) .

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

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

если с помощью какого-либо приближенного метода удается описать n+1

член разложения функции y y(x) в ряд Тейлора, то говорят, что метод n

-го порядка точности (в разложении в ряд Тейлора учитываются члены до порядка hn включительно). В методах Рунге-Кутта значение функции yi 1

в точке xi 1 определяется по формуле:

yi 1 yi yi

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

точности (или метод Эйлера), в котором для вычисления yi

2

используется касательная только в точке (xi , yi ) :

yi h f (xi ,yi ) .

Этот метод имеет довольно большую ошибку ограничения и часто оказывается неустойчивым - ошибка увеличивается с ростом x . В

исправленном методе Эйлера (метод Эйлера-Коши), который относится к семейству методов Рунге-Кутта второго порядка точности, вычисляется средний тангенс угла наклона касательной для двух точек (xi , yi ) и

(xi 1, yi h f (xi , yi )) :

y	(k1 k2 )	,k h f (x , y ),k			h f (x	, y k )
	2
i		1	i i	2	i 1	i	1 .

Одним из самых употребительных методов интегрирования дифференциальных уравнений является метод Рунге-Кутта четвертого порядка точности (или просто метод Рунге-Кутта без указания точности):

y

(k1 2 k2 2 k3 k4 )

,

i

6

h

k1

k h f (x , y ), k

2

h f (x

, y

),

1

i i

i

2

i

2

h

k2

k

3

h f (x

, y

), k

4

h f (x

, y

i

h k

).

i

2

i

2

i 1

3

Методы Рунге-Кутта переносятся на нормальные системы дифферен-циальных уравнений

yj '(x) f j (x, y1, y2 ,..., yn ), j 1,2,...,n ,

которые удобно записывать в векторной форме

Y '(x) F(x,Y ) ,

где Y (y1, y2 ,..., y n),F ( f1, f 2 ,..., f n) .

В расчетных формулах методов Рунге-Кутта y и f (x, y) заменяется

2

соответственно

(K1 (k11, k12,..., k1 n))

на Y и F(x,Y ) , коэффициенты k1 – на K1

, и т.д.

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

Например, из уравнения y '' f (x, y, y ') с начальными условиями y(x0) y0 ,

y '(x0) p0 и заменой	y1 y ,y2 y '	получим систему
	y	' y	,
	1	2
	y2 ' f (x, y1, y2 );

с начальными условиями y1 (x0 ) y0 , y2 (x0 ) p0 или в векторной форме:

																							Y ' F (x,Y ),Y (x 0) Y 0 ,
где Y (y1, y 2) , F ( f1, f 2 ) (y2 , f ) ,																											f f (x, y1, y 2) , Y0 (y0 , p0 ) .
		Формулы метода Рунге-Кутта для системы дифференциальных
уравнений второго порядка имеют вид ( в формулах индекс i																																								при y1i , y2i
опущен).
		В общем случае.																									Для уравнения 2-го порядка,
																											сведенного к системе 2-х
уравнений.
k11 h f1 (x i ,y1 ,y2 )																											h y2
k12 h f2 (x i ,y1 ,y2 )																											h f (xi ,y1,y 2)
k	21	h f (x				h				, y		k11			, y			2				k12			)		h ( y	2			k12			)
		1		i	2					1		2								2												2
k		h f		(x			h			, y			k11			, y							k12			)	h f (x								h	, y			k11		, y			k12	)
	22		2	i	2					1		2							2	2												i		2				1		2		2	2
k		h f (x			h			, y			k21			, y							k22			)			h ( y						k22				)
					2
	31	1		i						1		2					2			2										2				2

2