Квашнин С.Е. - Сборник научных программ на Паскале (1040523), страница 5

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

Текст из файла (страница 5)

Система линейных уравнений с ленточной матрицей [F4]Процедура bandsolve (band- лента, solve- решать) находит решение матричногоуравнения a*x=b, где а- матрица большого порядка, все нулевые элементыкоторой содержатся в узкой полоске (ленте), сконцентрированной вдоль главнойдиагонали. Параметр n- порядок матрицы а, m- ширина ленты (обязательнонечетное число).Процедура bandsolve весьма эффективна, потому что оперирует только сленточной частью матрицы а, заданной в массиве с[1:n,1:m]. Ленточныеэлементы i-й строки матрицы а размещаются в i-й строке матрицы с так, чтоэлемент a[i,j] переходит в элемент c[i,j-i+(m+1)/2].Все ленточные элементы матрицы а, как нулевые, так и ненулевые, должныприсудствовать в массиве с.

Значения неопределенных (указанным способом)элементов массива с, таких как c[1,1] или c[n,m], могут быть произвольными.Массив v вначале содержит вектор b. После выполнения процедуры bandsolve вмассиве v[1:n] содержится решение (вектор x). Значение массива с зати рается впроцессе решения, выполняемого методом исключения по Гауссу с перестоновкойстрок, после которой следует обратная подстановка.{ Algorithm 195b, adopted by S.Kwashnin, RL7, MSTU,Moscow, USSR28.11.89, 28.05.91}procedure bandsolve(var c1,v1 ; n,m : integer; eps : real;var signal: boolean);label sign;vart: realType;i,j,jm,lr,piv,r: integer;35c : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute c1;v : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute v1;Function Indx(n,i,j : integer):integer;begin Indx:=n*pred(i)+j end;beginsignal:=true;lr:=(m+1) div 2;for r:=1 to lr-1 dofor i:=1 to lr-1 dobeginfor j:=2 to m do c[Indx(n,r,j-1)]:=c[Indx(n,r,j)];c[Indx(n,n+1-r,m+1-i)]:=0; c[Indx(n,r,m)]:=0;end;for i:=1 to n-1 dobegin piv:=i;for r:=i+1 to lr doif abs(c[Indx(n,r,1)])>abs(c[Indx(n,piv,1)]) then piv:=r;if abs(c[Indx(n,piv,1)])< eps then begin signal:=false; goto sign end;if piv <> i thenbegin t:=v[i];v[i]:=v[piv]; v[piv]:=t;for j:=1 to m dobegin t:=c[Indx(n,i,j)];c[Indx(n,i,j)]:=c[Indx(n,piv,j)]; c[Indx(n,piv,j)]:=t;end;{j}end;{swi line}v[i]:=v[i]/c[Indx(n,i,1)];for j:=2 to m do c[Indx(n,i,j)]:=c[Indx(n,i,j)]/c[Indx(n,i,1)];for r:=i+1 to lr dobegin t:=c[Indx(n,r,1)];v[r]:=v[r]-t*v[i];for j:=2 to m do c[Indx(n,r,j-1)]:=c[Indx(n,r,j)]-t*c[Indx(n,i,j)];c[Indx(n,r,m)]:=0;end; {r}if lr<>n then lr:=lr+1end; {triangulation}v[n]:=v[n]/c[Indx(n,n,1)];jm:=2;for r:=n-1 downto 1 dobeginfor j:=2 to jm dov[r]:=v[r]-c[Indx(n,r,j)]*v[r-1+j];if jm<> m then jm:=jm+1;end;sign:end;{BandSolve195}Пример программы:program testBandSolve;Uses Crt;36const n=5; m=3;typemass_N = array[1..n] of real;mass_NM = array[1..n,1..m] of real;vareps: real;signal: boolean;i,j: integer;c: mass_NM;v: mass_N;{$I BANDSOLV.PAS}begineps:=1E-13;c[1,1]:=0; c[1,2]:=1; c[1,3]:=2;c[2,1]:=3; c[2,2]:=4; c[2,3]:=5;c[3,1]:=7; c[3,2]:=8; c[3,3]:=9;c[4,1]:=2; c[4,2]:=6; c[4,3]:=0;c[5,1]:=4; c[5,2]:=5; c[5,3]:=0;{ BandSolve(c,n,m,eps,signal,v);}BandSolve(c,v,n,m,eps,signal);writeln('signal=',signal);for i:=1 to n dowrite(v[i]);writeln;end.ПОДПРОГРАММА ПЕРЕМНОЖЕНИЯ КВАДРАТНОЙ МАТРИЦЫ НА ВЕКТОР(MultyMVm ) .MultyMV – эта подпрограмма предназначена для перемножения квадратнойматрицы Matr , размерность которой (n × n) ,на вектор столбец Vectразмерности (n×1) ,в результате чего вычисляется вектор столбец Rezразмерности (n×1) .MultyMV(n,nmax,Matr,Vect,Rez );Вводимые параметры:n – размерность перемножаемой матрицы,nmax – размерность матрицы в ее описании,Matr – квадратная матрица размерности (n × n) ;Vect – вводимый вектор столбец размерности (n×1);Rez – вектор столбец результат размерности (n ×1) .Программа используется совместно с модулем Kernel ,в котором приводитсяописание Realtype.37Текст программы (turbo pas)Procedure MultyMV( n:integer; var Matr1,Vect1,Rez1);var i,k:integer;Matr :array[1..(2*MaxInt)div sizeof(Realtype)] of Realtype absolute Matr1;Vect :array[1..(2*MaxInt)div sizeof(Realtype)] of Realtype absolute Vect1;Rez: array[1..(1*MaxInt)div sizeof(Realtype)] of Realtype absolute Rez1;Function Indx(n,i,j:integer):integer;beginIndx:=n*pred(i)+j;end;beginfor i:=1 to n dobeginRez[i]:=0;for k:=1 to nRez[i]:=Rez[i]+Matr[Indx(n,i,k)] *Vect[k];end;end ;Пример программыUses Kernel,MultyMVmU;const nmax=14;n=3;var i,j :integer;Y,yy:array[1..nmax] of RealType;f :array[1..nmax,1..nmax] of RealType;beginfor i:=1 to n dofor j:=1 to n dof[i,j]:=0;for i:=1 to n dof[i,j]:=3;for i:=1 to n do Y[i]:=i;MultyMVMax(n,nmax,f,Y,yy);for I:=1 to n do writeln(Y[i]:10:2,’ ‘,yy[i]:10:2);writeln;end.Описание к примеру:В этом примере квадратная матрица f, описанная какf : array [ 1..nmax , 1..nmax ] of realtype ,где nmax – ee размерноcть,умножается на вектор столбец Y,описанный какY : array [1:nmax] ,результат заносится в вектор столбец yy ,описанный так же как и вектор Y.Вызов подпрограммы осуществляется оператором MultyMVMax(n,nmax,f,Y,yy),в котором n указывает на фактический размер минора матрицы А размерности38[1:n,1:n] и подвектор V.Входной вектор 123Результат 369Ортогонализация промежуточных значений векторов решений системдифференциальных уравнений (ORTO){---------------------------13.04.91 -----------------------------------}{ КVASHNIN S.E.

}procedure ORTOS(var Inp1,Outp1,Omega1;n : integer; { размерность вектора }m : integer; { количество векторов однородной системы,m - всегда четное и равно n или n/2 }gomogen : boolean { true - если решается однородная система,false - если решается неоднородная система } );vari,j,k,kk,NumVect : integer;a: realType;Inp : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute Inp1;Outp : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute Outp1;Omega: array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute Omega1;Function Indx(n,i,j : integer):integer;{ n - кол-во столбцов }beginIndx:=n*pred(i)+j;end;function Scalar(var X1,Y1; ifirst, jsecond : integer): RealType;var k : integer;s : RealType;X : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute X1;Y : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute Y1;begins:=0;for k:=1 to n dos:=s+X[indx(n,ifirst,k)]*Y[indx(n,jsecond,k)];Scalar:=send;begin{ n, m = nh or n gomogen }if ( odd(n) or not( ( m <> n) xor (m <> n div 2)) ) thenbegin writeln(#7'Not correct call ORTOS !!!');Halt end;if gomogen then NumVect:=m else NumVect:=m+1;for i:=1 to NumVect dofor j:=1 to NumVect do Omega[indx(NumVect,i,j)]:=0;Omega[indx(NumVect,1,1)]:=sqrt(Scalar(Inp,Inp,1,1));for i:=1 to n doOutp[indx(n,1,i)]:=Inp[indx(n,1,i)]/Omega[indx(NumVect,1,1)];for k:=2 to m do39beginfor j:=1 to k-1 doOmega[indx(NumVect,j,k)]:=Scalar(Inp,Outp,k,j);a:=Scalar(Inp,inp,k,k);for j:=1 to k-1 do a:=a-sqr(Omega[indx(NumVect,j,k)]);if a < 0 then begin writeln(#7,'Error in ORTO '); halt end;Omega[indx(NumVect,k,k)]:=sqrt(a);for i:=1 to n dobeginOutp[indx(n,k,i)]:=Inp[indx(n,k,i)];for j:=1 to k-1 doOutp[indx(n,k,i)]:=Outp[indx(n,k,i)]Omega[indx(NumVect,j,k)]*Outp[indx(n,j,i)];Outp[indx(n,k,i)]:=Outp[indx(n,k,i)]/Omega[indx(NumVect,k,k)];end; { i }end; { k }if gomogen then exit;for j:=1 to m do{ ?? } Omega[indx(NumVect,j,NumVect)]:=Scalar(Inp,Outp,NumVect,j);for i:=1 to n dobeginOutp[indx(n,NumVect,i)]:=Inp[indx(n,NumVect,i)];for j:=1 to m doOutp[indx(n,NumVect,i)]:=Outp[indx(n,NumVect,i)]Omega[indx(NumVect,j,NumVect)]*Outp[indx(n,j,i)];end;Omega[indx(NumVect,NumVect,NumVect)]:=1;end; { ORTOS }{-----------------------------------------------------------------------------}program testOrto; { 10.05.90 Eleonora ?? , 15.04.91 Kvashnin S.

}uses Crt,FrazaXYU;const n=6; m=n div 2;Nmax=9;type RealType = double;mas_m_n = array[1..Nmax,1..Nmax] of RealType;mass_mm = array[1..Nmax,1..Nmax] of RealType;varOutpS,InpS,InpSC : mas_m_n;OmegaS: mass_mm;i,j,k: integer;Rez: array[1..n] of RealType;reler : double;ch1 : char;{$I ORTOSM.PAS}{ I MULTYMVM.PAS}beginClrScr;Randomize;Write('Test');40FrazaXY(7,1,'*Automatic or *Manual generated input Vectors [*A/*M] ?',ch1);WriteLn(' Array Inp:');for i:=1 to m dobeginGoToXY(3,i+2);Write(i,' : ');for j:=1 to n dobeginGoToXY(j*8+5,i+2);if ch1='M' then Read(InpS[i,j]) elsebeginInpS[i,j]:=random(1000);write(InpS[i,j]:8:2);end;end;end;writeln;ORtoSM(InpS,OutpS,OmegaS,Nmax,n,m+1,false);WriteLn(' Outp : ');for i:=1 to m dobeginWrite(' N = ',i,' ');for j:=1 to n doWrite(OutpS[i,j]:8:2);writeln;end;writeln(' Omega ');for i:=1 to m dobeginfor j:=1 to m doWrite(OmegaS[i,j]:8:2);writeln;end;reler:=0;for j:=1 to n dofor i:=1 to m dobeginInpSC[i,j]:=0;for k:=1 to m doInpSC[i,j]:=InpSC[i,j]+OmegaS[k,i]*OutpS[k,j];RELER:=reler+sqr(InpS[i,j]-InpSC[i,j]);end;WriteLn(' Input Test (must equal Inpt) : ');for i:=1 to m dobeginWrite(' N = ',i,' ');for j:=1 to n doWrite(InpSC[i,j]:8:2);writeln;41end;writeln('Относительная ошибка ортогонализации всех введенных');writeln(' векторов - RELER = ',reler:15:10);readln;end.Литература1.

Характеристики

Список файлов книги