Квашнин С.Е. - Сборник научных программ на Паскале (1040523), страница 2
Текст из файла (страница 2)
{ program }АЛГОРИТМ № 41б . Вычисление определителя [F3]Эта процедура вычисляет определитель матрицы а размерности n*n методомтриангуляции. В процессе вычисления матрица а не сохраняет свое значение.11Function Det(var a1; n: integer): RealType;{ From "CACM" Algorithm 41 b, Adapted by S.Kvashnin, August 1987 }label fin;vart,d,max : RealType;i,j,k : integer;a: array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute a1;Function Indx(n,i,j : integer):integer;begin Indx:=n*pred(i)+j end;begind:=1.0;for k:=1 to n dobegin max:=0;for i:=k to n dobegin t:=a[Indx(n,i,k)];if abs(t) > abs(max) thenbegin max:=t; j:=i end;end;{i}if max=0 then begin d:=0; goto fin end;if j<> k thenbegin d:=-d;for i:=k to n dobegin t:=a[Indx(n,j,i)]; a[Indx(n,j,i)]:=a[Indx(n,k,i)];a[Indx(n,k,i)]:=t;end;end;for i:=k+1 to n dobegin t:=a[Indx(n,i,k)]/max;for j:=k+1 to n doa[Indx(n,i,j)]:=a[Indx(n,i,j)]-t*a[Indx(n,k,j)]end; {i}d:=d*a[Indx(n,k,k)]end; {k}fin: det:=dend; {det}Пример программы :program Determinant;Const nd=3;var d : RealType;a : array[1..nd,1..nd] of RealType;{$I DET.PAS}begin12a[1,1]:=1.; a[1,2]:=-3; a[1,3]:=2.; a[1,4]:=0 ;a[2,1]:=-3.;a[2,2]:=3.; a[2,3]:=-1.; a[2,4]:=1.;a[3,1]:=2.; a[3,2]:=-1.;a[3,3]:=0.; a[3,4]:=-2.;a[4,1]:=1.; a[4,2]:=-3.;a[4,3]:=2.; a[4,4]:=1.;d:=det(a,nd);writeln('Det=',d,' Must be = -1.00000000');end.АЛГОРИТМ № 42б.
Обращение матрицы [F1]Процедура invert обращает квадратную матрицу matr порядка n с помощьюряда элементарных операций над строками матрицы matr, расширенной путемдополнения ее единичной матрицей. Случай вырожденной матрицы указываетсязначением s=1 после выполнения процедуры.procedure invert(var matr2,matr11; n: integer;var s: integer);{ Algorithm 42b, "CACM", adapted by S.Kvashnin, 25.07.88}label test0,fin;var t: RealType;i,j,k,m: integer;a: array[1..20,1..40] of RealType; {1..n,1..2*n}matr : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute matr2;matr1 : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute matr11;Function Indx(n,i,j : integer):integer;begin Indx:=n*pred(i)+j end;beginm:=2*n; s:=0;for i:=1 to n dofor j:=1 to m doif j<=n then a[i,j]:=matr[Indx(n,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];13for 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[Indx(n,i,j)]:=a[i,j+n];s:=0;fin: END; {invert}Пример программы:program test_invert;Uses Kernel;const n=4;typemas_n_n=array[1..n,1..n] of RealType;var i,j,s : integer;h1,h: mas_n_n;{$I INVERT.PAS}Beginfor i:=1 to n dofor j:=1 to n doh[i,j]:=1/(i+j-1.0);invert(h,h1,n,s);for i:=1 to n dobeginfor j:=1 to n dowrite(h1[i,j]:11:4,'writelnend;writeln(' s=',s);end.');АЛГОРИТМ № 60б.
Вычисление интеграла по Ромбергу [D1]bПроцедура-функция rombint вычисляет значения интеграла∫ f ( x )dxcaпогрешностью порядка 2*k+2, где k>=0. При увеличении k на единицу времявычисления, грубо говоря, удваивается.function rombint(a,b :real; k :integer; swiint : char): real;vard,s,h: real;i,j,m,n : integer;t: array[1..21] of real;beginif k> 20 then writeln(' K < 20 !!!!! ');14d:=b-a;t[1]:=(FexRomb(a,swiint)+FexRomb(b,swiint))/2;n:=1;for i:=1 to k dobegins:=0; n:=2*n;h:=d/n;for j:=1 to n div 2 dos:=s+FexRomb(a+(2*j-1)*h,swiint);t[i+1]:=(2*s/n+t[i])/2;m:=1;for j:=i downto 1 dobeginm:=4*m;t[j]:=t[j+1]+(t[j+1]-t[j])/(m-1)end; {j}end; {i}rombint:=t[1]*dend; {rombint}Пример программы:program test_rombint;var int : real;function FexRomb(x: real; swiint : char): real;beginFexRomb:=sin(x);end;{$I ROMBINT.PAS}beginint:=rombint(0,pi,5,'a');writeln(int);end.АЛГОРИТМ № 77б.
Интерполяция, дифференцирование и интегрированиефункций. [ D1, D4, E1]Процедура difint может в зависимости от способа обращения к ней выполнятьинтерполяцию, дифференцирование, или интегрирование функций однойперемен ной ,которые на всем или на интересующей нас части интервалаадекватно описываются с помощью двупараболической кривой.Процедура основана на использовании интерполяционной схемы Лагранжа,приспособленной для осредненных парабол второго порядка. По данному методувычисляется производная от функций, численно определенной в точках 1,2,3 и 4,путем проведения одной параболы через точки 1,2,3 и второй параболы черезточки 2,3 и 4.Тогда производная в точке 2 равна среднему значению аналити ческих производных двух парабол ,т.е.
коэффициенты параболы (a1*x2^2+b1*x2+c1),проходящей через точки 1,2 и 3, и коэффициенты параболы (a2*x2^2+b2*x2+c2),проходящей через точки 2,3 и 4, определяются с помощью уравнений Лагранжа,как показано ниже .Арифметические средние значения этих коэффициентов15a=(a1+a2)/2, b=(b1+b2)/2, c=(c1+c2)/2используются для получения производной на интервале от точки 2 до точки 3 поформуле (2*a*x+b).Интерполяция производится аналогично, за исключением того, что конечнойформулой является парабола (a*x^2+b*x+c). Интегрирование производитсяподобно построению кривой по точкам.
Например, интеграл между некоторымидвумя точками, скажем, 2 и 3,является осредненным интегралом двух параболмежду границами независимых координат для точек 2 и 3. Осреднениепроизводится для каждого интервала вдоль абсциссы, так как полученныерезультаты накапливаются для вычисления определенного интеграла.С помощью уравнений Лагранжа находятся коэффициенты a[k],b[k],c[k] путемвычисления t =yj /m∏ ( x j − x i ), y j =i=1, i≠ jf ( x j ), j = 1,2,.., m.Тогдаmmmmma = ∑ t ;b = ∑ t ×x ;c = ∑ t ×x .∑∏ki kij kijj = 1, j ≠ ii =1i =1j = 1, i ≠ ji =1В данной процедуре эти формулы применяются для m=3, k=1,2.Значения формальных параметров поясняются ниже:x-массив табличных значений аргумента, имеющий размерность [1:n];всезначения аргумента различны и образуют монотонно возрастающую последовательность;y-массив табличных значений функций, имеющий размерность [1:n];jt-параметр, задающий режим выполнения процедуры difint.
Для интерполяциинужно положить jt=1; для дифференцирования jt=2; а для интегрирования jt=3;xarg-значение аргумента, при котором нужно производить интерполяцию илидифференцирование. Если xarg находится вне интервала (x[1],x[n-1]), тоформулы Лагранжа работают как экстраполяционные. При интегрированиизначение xarg несущественно;x1,x2- нижний и верхний пределы интегрирования. При интерполяции идифференцировании значения x1 и x2 несущественны.program difint777;type mass_n77 = array[1..100] of real;varn,jt,i: integer;xarg,x1,x2,res,pi : real;x,y: mass_n77;procedure difint77(n,jt : integer; xarg,x1,x2 : real; var x11,y1;var res : real);label start,no5,no6,no9,intrp,diff,intgr,no17,term,lagr,fin;varca,cb,cc,a,b,c,s1,s2,t1,t2,t3,sum : real;j,js,j2,i,j1,jj2: integer;y : array[1..(2*MaxInt) div SizeOf(double)] of double absolute y1;x : array[1..(2*MaxInt) div SizeOf(double)] of double absolute x11;beginstart: if jt=3 then goto intgr;16if xarg>=x[n-1] thenbegin j:=n-1; js:=1; goto term end;if xarg<=x[2] thenbegin j:=2; js:=1; goto term end;{}js:=2;for i:=2 to n dobegin if xarg<x[i] then goto term; j:=i end;no5: ca:=a; cb:=b; cc:=c;js:=3; j:=j+1;goto term;no6: a:=(ca+a)/2; b:=(cb+b)/2; c:=(cc+c)/2;no9: if jt=2 then goto diff;intrp: res:=a*sqr(xarg)+b*xarg+c;goto fin;diff: res:=2*a*xarg+b;goto fin;intgr: sum:=0; s1:=x1;j2:=n; j1:=2;for i:=1 to n dobeginif x1<=x[i] then goto no17;j1:=j1+1;end;no17: for i:=1 to n dobegin j2:=j2-1;if x2>=x[j2+1] then goto lagr;end;term: j1:=j; j2:=j;lagr: if j2 <= j1 then jj2:=j1+1 else jj2:=j2;for j:=j1 to jj2 dobegin a:=x[j-1]-x[j];b:=x[j-1]-x[j+1]; c:=a-b;t1:=y[j-1]/(a*b);t2:=y[j]/(a*c);t3:=-y[j+1]/(b*c);a:=t1+t2+t3;b:=-(x[j]+x[j+1])*t1-(x[j-1]+x[j+1])*t2-(x[j-1]+x[j])*t3;c:=x[j]*x[j+1]*t1+x[j-1]*x[j+1]*t2+x[j-1]*x[j]*t3;if jt<> 3 thenbeginif js=1 then goto no9;if js=2 then goto no5;if js=3 then goto no6;end;if j=j1 thenbegin ca:=a; cb:=b; cc:=c; end elsebeginca:=(a+ca)/2; cb:=(b+cb)/2; cc:=(c+cc)/2;end;s2:=x[j];sum:=sum+ca*(s2*s2*s2-s1*s1*s1)/3+cb*(s2*s2-s1*s1)/2+cc*(s2-s1);ca:=a; cb:=b; cc:=c;17s1:=s2end; {j}res:=sum+ca*(x2*x2*x2-s1*s1*s1)/3+cb*(x2*x2-s1*s1)/2+cc*(x2-s1);fin: end; {difint77}beginpi:=3.141593;write('n=');readln(n);for i:=0 to n dobeginx[i+1]:=pi*i/n;y[i+1]:=sin(i*pi/n);end;jt:=3;x1:=0;x2:=pi;difint77(n,jt,xarg,x1,x2,x,y,res);writeln('Res=',res,' Must be = 2.0000000');write('n,jt,xarg= '); readln(n,jt,xarg);x1:=0; x2:=2;for i:=1 to n+1 dobeginx[i]:=x1+(i-1)*(x2-x1)/n;y[i]:=exp(x[i]);end;difint77(n,jt,xarg,x1,x2,x,y,res);writeln('Res=',res);end.АЛГОРИТМ №120б.
Обращение матрицыЭтот алгоритм выполняет обращение матрицы а с записью результата на местоматрицы а, имеющей размерность [1:n,1:n]. Если в процессе вычислениянекоторый главный элемент имеет абсолютное значение меньше eps, топроисходит переход к метке signal. Переменная det будет равна значениюопределителя исходной матрицы при нормальном выходе из процедуры и равнанулю или очень малому числу при выходе к метке signal.{--------------------Adapted by Kvashnin S.E. , Pavlova ----------- }{-----------------------26.12.89------------------------------------}procedure inversion2(n : integer; eps : real; var signal : boolean;var a1; var det : RealType);label rpt,nexti;vari,j,k,p,r : integer;y,w: RealType;b,c : array[1..40] of RealType; { array[1..n] }z: array[1..40] of integer;a: array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute a1;Function Indx(n,i,j : integer):integer;begin Indx:=n*pred(i)+j end;18beginsignal:=false;det:=1.0;for j:=1 to n do z[j]:=j;for i:=1 to n dobegin k:=i; y:=a[indx(n,i,i)];r:=i-1; p:=i+1;for j:=p to n dobegin w:=a[indx(n,i,j)];if abs(w) > abs(y) thenbegin k:=j; y:=w endend;{ j }det:=y*det;if abs(y) < eps then exit;y:=1.0/y;for j:=1 to n dobegin c[j]:=a[indx(n,j,k)]; a[indx(n,j,k)]:=a[indx(n,j,i)];a[indx(n,j,i)]:=-c[j]*y;a[indx(n,i,j)]:=a[indx(n,i,j)]*y; b[j]:=a[indx(n,i,j)];end; { j }a[indx(n,i,i)]:=y;j:=z[i]; z[i]:=z[k]; z[k]:=j;for k:=1 to n dofor j:=1 to n doif not ((k=i) or (j=i)) then a[indx(n,k,j)]:=a[indx(n,k,j)]-b[j]*c[k]end;{ i;}for i:=1 to n dorpt: begin k:=z[i];if k=i then goto nexti;for j:=1 to n dobegin w:=a[indx(n,i,j)]; a[indx(n,i,j)]:=a[indx(n,k,j)];a[indx(n,k,j)]:=wend; { j }p:=z[i]; z[i]:=z[k];z[k]:=p;det:=-det;goto rpt;nexti: end; { i }signal:=true;end; { inversion2 }АЛГОРИТМ №126б.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
