Пример программы (лабораторная работа 3)

Варианты задания.

Исходная СЛАУ

,

где – номер факультета; – номер группы; – номер студента по журналу.

Не выполняется условие диагонального преобладания для 1-го и 2-го уранений. Их перестановка позволяет выполнить это условие:

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

– норма по модулю.

1. Обычная формулировка

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

integer,parameter::n=3,kmax=100 real,parameter::s=2,g=3,r=s-g,k=1,r1=s+g+k,a33=2*r1,eps=0.001 real,dimension(n,n)::A=(/4.,1.,s,-1.,5.,g,1.,1.,a33/) real,dimension(n)::b, xk,xk1 b(1)=r+4*k-2;b(2)=r+k+10;b(3)=2*(r*r1+g)+k*s print 100,'matrix A',((a(i,j),j=1,n),i=1,n);print 100,'vector b',b 100 format(/10x,a/(3f10.1)) print 100,'SIMPLE ITERATION'; xk=0; do ki=0,kmax; zk=0; xk1=b; do i=1,n; do j=1,n; if(i/=j) xk1(i)=xk1(i)-a(i,j)*xk(j); enddo; xk1(i)=xk1(i)/a(i,i);zk=zk+abs(xk1(i)-xk(i)); enddo; xk=xk1; if(zk<eps)goto 3 enddo; 3 print '(2x,"ki=",i3,3x,"zk=",f8.6)',ki,zk; print 100,'result',xk print 100,'SEIDEL'; xk=0; do ki=0,kmax; zk=0; do i=1,n; y=b(i) do j=1,n; if(i/=j) y=y-a(i,j)*xk(j); enddo; y=y/a(i,i);zk=zk+abs(y-xk(i)); xk(i)=y enddo; if(zk<eps)goto 4 enddo 4 print '(2x,"ki=",i3,3x,"zk=",f8.6)',ki,zk; print 100,'result',xk   end

Результаты счета

matrix A

4.0 -1.0 1.0

1.0 5.0 1.0

2.0 3.0 12.0

vector b

1.0 10.0 -4.0

SIMPLE ITERATION

ki= 6 zk=0.000065

result

1.0 2.0 -1.0

SEIDEL

ki= 4 zk=0.000827

result

1.0 2.0 -1.0

Press any key to continue

2. Maтричная формулировка

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

use dfimsl use linear_operators integer,parameter::n=3,kmax=100 real,parameter::s=2,g=3,r=s-g,k=1,r1=s+g+k,a33=2*r1,eps=0.001 real,dimension(n,n)::A=(/4.,1.,s,-1.,5.,g,1.,1.,a33/),L,U,C1,C2,D real,dimension(n)::b,b1,b2,xk,xk1,di b(1)=r+4*k-2;b(2)=r+k+10;b(3)=2*(r*r1+g)+k*s print 100,'matrix A',((a(i,j),j=1,n),i=1,n);print 100,'vector b',b 100 format(/10x,a/(3f10.1))   D=diag(diagonals(A)); L=0;U=0;do i=1,n-1;do j=i+1,n;U(i,j)=A(i,j);enddo;enddo;L=A-D-U; print 100,'matrix L',((L(i,j),j=1,n),i=1,n); print 100,'matrix D',((d(i,j),j=1,n),i=1,n); print 100,'matrix U',((U(i,j),j=1,n),i=1,n); xk=0 print 100,'SIMPLE ITERATION'; C1=D.ix.(A-D); b1=D.ix.b do ki=0,kmax xk1=b1-matmul(C1,xk) zk=sum(abs(xk1-xk)) xk=xk1;if(zk<eps)goto 1 enddo 1 print '(2x,"ki=",i3,3x,"zk=",f8.6)',ki,zk print 100,'result',xk1 xk=0 print 100,'SEIDEL'; C2=(L+D).ix.U;b2=(L+D).ix.b do ki=0,kmax xk1=b2-matmul(C2,xk) zk=sum(abs(xk1-xk)) xk=xk1;if(zk<eps)goto 2 enddo 2 print '(2x,"ki=",i3,3x,"zk=",f8.6)',ki,zk print 100,'result',xk1   end

Замечание. Здесь приняты следующие обозначения

,

где

– нижняя треугольная матрица,

– диагональная матрица,

– верхняя треугольная матрица.

Простая итерация: ,

где .

Метод Зейделя: ,

где ,

Операторы

use dfimsl

use linear_operators

означают подключение библиотек стандартных подпрограмм dfimsl и
linear_operators, соответственно.

В библиотеке linear_operators возможно использовать следующие матричные операции

Математика	Fortran 90
	A.x.B
	.i.A
	A.ix.B
	A.xi.B

Результаты счета

matrix A

4.0 -1.0 1.0

1.0 5.0 1.0

2.0 3.0 12.0

vector b

1.0 10.0 -4.0

matrix L

0.0 0.0 0.0

1.0 0.0 0.0

2.0 3.0 0.0

matrix D

4.0 0.0 0.0

0.0 5.0 0.0

0.0 0.0 12.0

matrix U

0.0 -1.0 1.0

0.0 0.0 1.0

0.0 0.0 0.0

SIMPLE ITERATION

ki= 6 zk=0.000065

result

1.0 2.0 -1.0

SEIDEL

ki= 4 zk=0.000827

result

1.0 2.0 -1.0

Press any key to continue

Дата добавления: 2014-01-11