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

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

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

Решение системы линейных уравнений методом ГауссаЭта процедура предназначена для решения системы u линейныхалгебраических уравнений методом последовательного исключения неизвестных.a[1:u,1:u+1]- расширенная матрица системы ;u- число неизвестных;y[1:u]- вектор решения;Если система не имеет решения или имеет множество решений тоосуществляется выход из процедуры к метке signal126.19Procedure Gauss126(u:integer; var a1,y1; var signal:boolean);label 10,11,12;var temp : RealType;i,j,k,m,n: integer;a : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute a1;y : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute y1;Function Indx(n,i,j : integer):integer;{ n - кол-во столбцов }begin Indx:=n*pred(i)+j end;beginn:=0; signal:=true;10: n:=n+1;for k:=n to u doif a[Indx(u+1,k,n)]<>0 then goto 11;signal:=false; exit;11: if k=n then goto 12;j:=u+1;for m:=n to j dobegin temp:=a[Indx(u+1,n,m)]; a[Indx(u+1,n,m)]:=a[Indx(u+1,k,m)];a[Indx(u+1,k,m)]:=temp end;12: for j:=u+1 downto n doa[Indx(u+1,n,j)]:=a[Indx(u+1,n,j)]/a[Indx(u+1,n,n)];m:=u+1;for i:=k+1 to u dofor j:=n+1 to m doa[Indx(u+1,i,j)]:=a[Indx(u+1,i,j)]-a[Indx(u+1,i,n)]*a[Indx(u+1,n,j)];if n<>u then goto 10;for i:=u downto 1 dobegin y[i]:=a[Indx(u+1,i,m)];for k:=i-1 downto 1 doa[Indx(u+1,k,m)]:=a[Indx(u+1,k,m)]-a[Indx(u+1,k,i)]*y[i]end;{ i }END;Пример программы :program Gauss126;{ this ia a program of Gauss method's for solving the sistem of algebr.equations , see book CACM (USA), Algorithm numb.126 b }label 1000,1126;const uu=4;uuu=uu+1;type mass=array[1..uu] of double;mass2D=array[1..uu,1..uuu] of double;realType = double;var a: mass2D;yy: mass;uv: integer;sig: boolean;20{$I GAUSS.PAS}{ little example}begina[1,1]:=1; a[1,2]:=0.17; a[1,3]:=-0.25; a[1,4]:=0.54;a[1,5]:=0.3;a[2,1]:=0.47; a[2,2]:=1; a[2,3]:=0.67; a[2,4]:=-0.32;a[2,5]:=0.5;a[3,1]:=-0.11; a[3,2]:=0.35; a[3,3]:=1;a[3,4]:=-0.74;a[3,5]:=0.7;a[4,1]:=0.55;a[4,2]:=0.43;a[4,3]:=0.36;a[4,4]:=1;a[4,5]:=0.9;uv:=4;Gauss126(uv,a,yy,sig);writeln('', yy[1]:9:6,' ',yy[2]:9:6,' ',yy[3]:9:6,' ',yy[4]:9:6,' ',sig);writeln('Must be - 0.440888 -0.363031 1.166798 0.393567 ');1000: end.АЛГОРИТМ №163б.

Модифицированная функция Ханкеля [S17]Процедура expk (exponential- экспоненциальная, а k соответствует Kp(x))вычисляет модифицированную функцию Ханкеля exp(x)*Kp(x) с заданнойточностью e по интегральной формуле∞e K p = ∫ e x (1− cht ) × ch( pt ) dtx0program expk;uses Crt;function stp(z,p:real):real;var a,b:real;beginb:=p*ln(z);stp:=exp(b);end;label1,2;varexk,f,g,h,r,s,u,z,zp,p,x,e :double;beginClrScr;read(p,x,e);r:=0; h:=1;1:g:=r; s:=0;z:=exp(0.5*h); u:=z*z;2:zp:=stp(z,p);21f:=0.5*exp(x*(1-0.5*(z+1/z)))*(zp+1/zp);s:=s+f; z:=z*u;if f >= e then goto 2;r:=h*s; h:=0.5*h;if abs(r-g) >= e then goto 1;exk:=r;window(10,5,70,20);textbackground(15);textcolor(0);ClrScr;writeln('REZULT ');writeln(^G,' K_p:=',(2/PI*exk):5:5);end.Пример программы:program hankel;const pi=3.14159;var a,p1,x1,e1,ccc :real;m: integer;{$I 163pas.pas}beginwriteln;read (x1);p1:=0;e1:=0.0001;a:=2*expk(p1,x1,e1)/pi/exp(x1);writeln(a);readln;END.АЛГОРИТМ № 170б.

Определитель с полиномиальными элементами [F3]Процедура polimatrix раскрывает определитель общего вида, в котором каждыйэлемент является полиномом.Эта программа полезна для исследования задачдинамической устойчивости при использовании аппроксимации преобразующейфункции.Здесь сначала выполняется один из процессов триангуляризацииполиномиальной матрицы с вещественными коэффициентами, затем умножениемдиагональных элементов формируется детерминантный полином. Указаннаяздесь полиномиальная матрица имеет в качестве элементов полиномы видаn∑ ai x ii=0После триангуляризации все элементы ниже главной диагонали равнынулю.Далее при раскрытии определителя там формируются ненулевые члены.

Врезультате можно проверять, например, критерий устойчивости путемвычисления корней сформированного таким образом характеристическогоуравнения, пользуясь какой-нибудь подходящей программой нахождения корней.22Рассмотрим для примера полиномиальную матрицу с квадратичнымиэлементами (n=2). В этом случае массив а должен иметь размерность[1:p,1:p,1:m], где p- порядок матрицы, а m=n*p+1. Здесь первый индекссоответствует строке, второй- столбцу, а третий- коэффициенту полинома.Следовательно, перед обращением к процедуре постоянный член главногополиномиального элемента содержится в a[i,j,1], линейный- в a[i,j,2],аквадратичный- в a[i,j,3] и т.д.После выполнения программы коэффициенты детерминантного полиномасодержатся в массиве c[1:m].

Постоянный коэффициентбудет в c[1], линейный-вc[2], а квадратичный- в c[3] и т.д. переменная r будет равна числу коэффициентов детерминантного полинома. В общем случае r<>m, так как при раскрытииопределителя некоторые коэффициенты могут обратиться в нуль,то приобращении к процедуре нужно задавать r=m. Если полиномы, представ-ляющиесобой элементы матрицы, не все одного и того же порядка, то перед обращениемк процедуре нужно брать n равным наивысшей из степеней полиномов.Значение eps можно задавать, например, равным 2 10 = 8.В некоторых случаях неправильного хода решения осуществляется выход изпроцедуры к глобальной метке signal 170.{ 170b Kvashnin S.E., modify 26.05.91 }procedure polymatrix(var a1; p,n : integer; eps : RealType;{ rezults - }var r : integer; var c11; var signal : boolean);label n0,n2,n3,n4,n5,n6,n10,n12,n11,n14,n15;varsa,sb: RealType;i,j,k,j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,m : integer;c1,c2: array[1..21] of RealType;mat:array[1..10,1..10] of integer;a : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute a1;c : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute c11;Function Indx3(n,m,i,j,k : integer):integer;begin{ a[1..q,1..n,1..m] } { a[1..p,1..p,1..m] }Indx3:=n*m*pred(i)+m*pred(j)+k;end;beginsignal:=false;m:=n*p+1;for i:=1 to p dofor j:=1 to p dobegin j1:=0;for k:=1 to m doif a[Indx3(p,m,i,j,k)]<>0 then j1:=k;mat[i,j]:=j1end;{i}j1:=1;n0:j9:=0;for i:=j1 to p dobegin23if mat[i,j1]<0 then exit;if mat[i,j1]>0 thenbegin j9:=j9+1; j3:=i endend; {i}if j9<1 then exit;if j9>1 then goto n2;if j3<j1 then exit;if j3<> j1 thenfor j:=j1 to p dobegin j2:=mat[j3,j]; j4:=j2;if mat[j1,j]>j2 then j2:=mat[j1,j];mat[j3,j]:=mat[j1,j]; mat[j1,j]:=j4;for k:=1 to j2 dobegin sa:=a[Indx3(p,m,j3,j,k)];a[Indx3(p,m,j3,j,k)]:=a[Indx3(p,m,j1,j,k)];a[Indx3(p,m,j1,j,k)]:=saend;{k}end;{j}goto n12;n2:j3:=j1+1;for i:=j3 to p dobeginn3:if mat[i,j1]<0 then exit;if mat[i,j1]=0 then goto n11;if mat[j1,j1]<0 then exit;if (mat[j1,j1]<>0) and (mat[i,j1]>=mat[j1,j1]) then goto n5;n4:for j:=j1 to p dobeginj2:=mat[j1,j];j4:=j2;if mat[i,j]>j2 then j2:=mat[i,j];mat[j1,j]:=mat[i,j]; mat[i,j]:=j4;for k:=1 to j2 dobegin sa:=a[Indx3(p,m,i,j,k)];a[Indx3(p,m,i,j,k)]:=a[Indx3(p,m,j1,j,k)]; a[Indx3(p,m,j1,j,k)]:=-saend;end; {j}goto n3;n5:j7:=mat[i,j1]; j5:=mat[j1,j1];j6:=j7-j5; sb:=a[Indx3(p,m,i,j1,j7)]/a[Indx3(p,m,j1,j1,j5)];if abs(sb)<4 then goto n6;if j6<0 then exit;if j6=0 then goto n4;n6:for j:=j1 to p dobegin j5:=mat[j1,j];for k:=1 to j5 dobegin j7:=k+j6;if j7>m then goto n10;sa:=a[Indx3(p,m,i,j,j7)]-sb*a[Indx3(p,m,j1,j,k)];if abs(sa)<eps then a[Indx3(p,m,i,j,j7)]:=0else a[Indx3(p,m,i,j,j7)]:=saend; {k}end; {j}n10:for j:=j1 to p dobegin j7:=mat[j1,j]+j6;if mat[i,j]>=j7 then j7:=mat[i,j];24mat[i,j]:=0;for k:=1 to j7 doif a[Indx3(p,m,i,j,k)]<>0 then mat[i,j]:=kend; {j}n11:end; {i}goto n0;n12:j1:=j1+1;if j1<p then goto n0;for j:=1 to p dobegin j2:=mat[j,j];for k:=1 to j2 do c1[k]:=a[Indx3(p,m,j,j,k)];if j=1 then goto n14;for k:=1 to r do c2[k]:=c[k];for k:=1 to m do c[k]:=0;if j2<0 then exit;if j2=0 then goto n15;for k:=1 to j2 dofor j10:=1 to r dobegin j11:=k+j10-1;c[j11]:=c[j11]+c1[k]*c2[j10]end; {k}r:=j11;goto n15;n14:for k:=1 to j2 do c[k]:=c1[k];r:=j2;n15:end; {j}signal:=true;end; {polymatrix}Пример программы:program test_plmt;{ Algorithm 170b "CACM", adopted by S.Kvashnin, 25.07.88, 26.05.91n - max.

order of polynome into the i,j-element's matr. Ap - order matr. Am = n*p+1Befor begining: r:=m !!!a[i,j,1] - const.element of polynome,a[i,j,2] - linear element of polynome,a[i,j,3] - squear element of polynome etc.After ending: r - number elements of rezults polynom }Uses Crt;const p=2; n=4; m=9;type RealType = double;mas_3 = array[1..p,1..p,1..m] of RealType;mas_m = array[1..m] of RealType;var r,i,j,k : integer;a: mas_3;c: mas_m;sig: boolean;ch1 : string[1];25function sign(x: realType):integer;beginif x=0 then sign:=0 else sign:=round(x/abs(x));end;{$I POLYMATR.PAS}beginClrScr;writeln('Test PoliMatr 26.05.91 ') ;r:=m;for i:=1 to p dofor j:=1 to p dofor k:=1 to m doa[i,j,k]:=0;a[1,1,1]:=1;a[1,1,2]:=1;a[1,1,3]:=1;a[1,1,4]:=1;a[1,1,5]:=1;a[1,2,1]:=3;a[1,2,2]:=-2;a[1,2,4]:=-1;a[2,1,1]:=-1;a[2,1,2]:=1;a[2,2,1]:=1;polymatrix(a,p,n,1e-8,r,c,sig);gotoXY(1,7); write('Y(X)= ');for i:=1 to r dobeginif sign(c[i]) >=0 then ch1:='+' else ch1:='';write(ch1:1,c[i]:7:4,'*X^',i-1:1,' ');end;Write(' Must be y(x) = 4-4x+3x^2+2x^4 ');readln;end.АЛГОРИТМ №176 б.

Аппроксимация последовательности точек линейнойкомбинацией любых заданных функций [E2]Процедура surfit по данной последовательности m ординат и соответству ющихим значений от n наперед заданных основных функций f[i] одного или несколькихлинейно-независимых переменных аппроксимирует (в смысле методанаименьших квадратов) точки с помощью функции, имеющей формуa1*f1+a2*f2+...+an*fn, где ai- искомые коэффициенты.Вычисляются также вектор разностей ej и их средняя длина r.

Должны бытьпредусмотренны также веса, соответствующие данным точкам. Алгоритмсводится по существу к решению матричного уравнения f^t*w*f*a=f^t*w*z, гдеа- вектор искомых коэффициентов ai;w- вектор, содержащий диагональные элементыматрицы размерностью m*m из весов данных точек;z- вектор значений ординат;f- матрица размерностью m*n из соответствующихзначеий функции.Процедура surfit использует процедуру-параметр Invert, заменяющуювещественную матрицу gg (рабочий массив) на обратную.

Размерность массивовпараметров: a[1:n], e,w,z[1:m], f[1:m,1:n].26procedure invert(matr : mas_n_n; n: integer; var matr1 : mas_n_n;var s: integer);{ Algorithm 42b, "CACM", adopted by S.Kvashnin, 25.07.88}label test0,fin;var t: real;i,j,k,m: integer;a: array[1..20,1..40] of real; {1..n,1..2*n}beginm:=2*n; s:=0;for i:=1 to n dofor j:=1 to m doif j<=n then a[i,j]:=matr[i,j] elseif j=n+i then a[i,j]:=1 else a[i,j]:=0;for i:=1 to n dobegin k:=i;test0:if a[k,i]=0 thenbegin s:=1;if k<n then k:=k+1 else goto fin;goto test0end;if s=1 thenfor j:=1 to m dobegint:=a[k,j]; a[k,j]:=a[i,j]; a[i,j]:=tend;for j:=m downto i do a[i,j]:=a[i,j]/a[i,i];for k:=1 to n doif k<>i thenfor j:=m downto i doa[k,j]:=a[k,j]-a[i,j]*a[k,i]end; {i}for i:=1 to n dofor j:=1 to n do matr1[i,j]:=a[i,j+n];s:=0;fin: END; {invert}{ procedure invert(matr : mas_n_n; n: integer; var matr1 : mas_n_n;var s: integer);Algorithm 42b, "CACM", adapted by S.Kvashnin, 25.07.88}procedure surfit(f: msurMN; z,w: msurM; m,n : integer; var a:msurN;var e: msurM; var r: real);var gg, ggend : mas_n_n;i,j,k,ss : integer;s,g: real;beginfor i:=1 to n dofor j:=1 to n do27begins:=0;for k:=1 to m dos:=s+f[k,i]*f[k,j]*w[k];gg[i,j]:=send;invert(gg,n,ggend,ss);for i:=1 to n dofor j:=1 to n dogg[i,j]:=ggend[i,j];for i:=1 to n dobegins:=0;for j:=1 to m dobeging:=0;for k:=1 to n dog:=g+gg[i,k]*f[j,k];s:=s+g*z[j]*w[j]end;a[i]:=send; {i}s:=0;for i:=1 to m dobeging:=z[i];for j:=1 to n dog:=g-a[j]*f[i,j];s:=s+g*g; e[i]:=gend; {i}r:=sqrt(s/m);end; {sutfit}Пример программы:program surfitmain;Uses Crt;const n=2; m=6;typemsurMN = array[1..m,1..n] of real;msurM = array[1..m] of real;msurN = array[1..n] of real;mas_n_n = array[1..n,1..n] of real;varf: msurMN;z,w,e : msurM;a: msurN;i,j: integer;r,x: real;alfa,c,h,freq : real;{------------------------------------}{$I SURFIT.PAS}{$I KRILOV.PAS}{------------------------------- 13.10.89 ----------------------------------}28beginfor i:=1 to m do w[i]:=1.0;c:=5000; freq:=25000;alfa:=2*pi*freq/c;x:=0;h:=0.005;for i:=1 to m dobeginx:=x+h;z[i]:=2*cos(alfa*x);end;{ z[1]:=0.5; z[2]:=0.6; z[3]:=0.7;z[4]:=0.8; z[5]:=0.9; z[6]:=1.0;}x:=0.4;for i:=1 to m dobeginx:=x+0.1;f[i,1]:=sin(x);f[i,2]:=sh(x);end;surfit(f,z,w,m,n,a,e,r);ClrScr;writeln(a[1],a[2],r);writeln;for i:=1 to m do write(' ',e[i]);end.АЛГОРИТМ №178б.

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

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