Квашнин С.Е. - Сборник научных программ на Паскале

1МГТУ им. Н.Э. БауманаКафедра БМТ-2Сборникнаучных программ на ПаскалеУчебное пособие по курсовому проектированию для студентовфакультета БМТпод ред. Квашнина С.Е.Москва20052ВведениеВ учебном пособии представлены тексты научных программ из сборниковCASM [1-4]и переведенные с алгоритмического языка Алгол-60 на язык Паскаль, атакже ряд программ из книги Дж.Форсайта [5] и переведенные с языпа Фортран-4на Паскаль. Некоторые из программ написаны редактором настоящего сборника(перемножение матриц, процедура ортогонализации – ORTO).Все представленные в настоящем учебном пособии программы многократнотестировались и проверялись при решении разнообразных задач проектированиямедицинской техники.В настояший сборник вошли лишь наиболее надежные варианты программ,прошедшие испытания временем.Ко всем стандартным программам даны тестовые примеры на языке Паскаль, спомощью которых можно получить представление о работоспособностиалгоритма.

Нумерация алгоритмов дана по сборникам CASM. Часть идей, поадресации массивов заимствована из книги В.В.Фаронова [6].АЛГОРИТМ № 3б. Нахождение комплексных корней полинома методомБерстоу-Хичкока [C2]В процедуре bairstow итерация Берстоу-Хичкока используется дляпоследовательного нахождения пар корней полиномиального уравнения степениn с коэффициентами a[i] (i=0,1,..,n),где a[n]-свободный член.Выходные параметры:num-число пар найденных корней;x[i],y[i] (i=1,2,..,num)-пара действительных корней, если nat[i]=1,илидействительная и мнимая части комплексной пары, если nat[i]=-1;ex[i] указывает,по какому из следующих четырех условий произошелвыход из итерационального цикла при нахождении пары корней x[i],y[i];ex[i]=1,если u1,u1, стали по модулю меньше eps1;ex[i]=2,если поправки corp,corg становятся по модулю меньше eps2;ex[i]=3,если отношения corp/p,corg/g становятся меньше eps3;ex[i]=4,если чило итераций становится равным max.В последнем случае найденная пара корней недостоверна.При этом неделается никакой попытки найти остальные корни.Eps0 используется как нижняя граница знаменателя (det)в выражениях, изкоторых были найдены corp и corg.Входные массивы x,y,nat и ex должны иметь размерность [1:(n+1)/2].procedure bairstow(n:integer; var a1,x1,y1,nat1,ex1; eps0,eps1,eps2,eps3 :real;max: integer;var num : integer);label out,fin,lab1,lab2;var3p,q,corq,corp,u0,u1,v0,v1,w0,w1,det,s,t : RealType;i,j,k,m,n2,r,ii: integer;b,c{[1..n+1]}: array[0..20] of RealType;a : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute a1;x : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute x1;y : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute y1;nat: array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute nat1;ex : array[1..(2*MaxInt) div SizeOf(RealType)] of RealType absolute ex1;Function Indx(n,i,j : integer):integer;begin Indx:=n*pred(i)+j end;beginfor i:=0 to n do b[i]:=a[i];b[n+1]:=0; n2:=(n+1) div 2;r:=2*n2;for m:=1 to n2 dobeginp:=0; q:=0;for k:=0 to max dolab1:beginfor i:=0 to r do c[i]:=b[i];ii:=1;i:=r-2;repeatfor j:=0 to i dobegin c[j+1]:=c[j+1]-p*c[j];c[j+2]:=c[j+2]-q*c[j]end;i:=r-4; ii:=ii+1;until ii=3;u0:=c[r]; u1:=c[r-1];v0:=c[r-2]; v1:=c[r-3];w0:=-q*v1; w1:=v0-p*v1;det:=v0*w1-w0*v1;if abs(det)<eps0 thenbegin p:=p+1; q:=q+1;goto lab1end;corp:=(v0*u1-u0*v1)/det;corq:=(u0*w1-w0*u1)/det;p:=p+corp;q:=q+corq;if(abs(u0)<eps1) and (abs(u1)<eps1) thenbegin ex[m]:=1; goto lab2 end;if(abs(corp)<eps2) and (abs(corq)<eps2) thenbegin ex[m]:=2; goto lab2 end;if(abs(corp)<eps3*abs(p)) and (abs(corq)<eps3*abs(q))thenbegin ex[m]:=3; goto lab2 end;end; {k}ex[m]:=4;4lab2:s:=-p/2; t:=s*s-q;if t>=0 thenbegin nat[m]:=1; t:=sqrt(t);x[m]:=s+t; y[m]:=s-t;endelsebegin nat[m]:=-1; x[m]:=s;y[m]:=sqrt(-t)end;if ex[m]=4 then goto out;r:=r-2;for j:=0 to r dobegin b[j+1]:=b[j+1]-p*b[j];b[j+2]:=b[j+2]-q*b[j]end;if r<1 thenout:begin num:=m; goto fin end;if r<3 thenbegin m:=m+1; ex[m]:=1;p:=b[1]/b[0]; q:=b[2]/b[0];goto lab2;end;end; {m}fin:end; {bairstow}Пример программы:program test_brst;const n=6; num=3;type mas_0n=array[0..n] of real;imas_num=array[1..num] of integer;mas_num=array[1..num] of real;var i,j,numm: integer;a: mas_0n;x,y: mas_num;nat,ex: imas_num;e0,e1,e2,e3: real;{$I BAIRSTOW.PAS}Begine0:=1e-6; e1:=1e-6; e2:=1e-6; e3:=1e-6;a[0]:=1;a[1]:=-2;a[2]:=2;a[3]:=1;a[4]:=6;5a[5]:=-6;a[6]:=8;bairstow(n,a,e0,e1,e2,e3,10000,numm,x,y,nat,ex);for i:=1 to numm dowriteln(x[i],y[i],' ',nat[i],' ',ex[i]);End.АЛГОРИТМ № 4б.

Нахождение корней непрерывной функции методомделения интервала пополам [C5].Процедура bisec вычисляет функцию на концах вещественного интервала.Если нет перемены знака функции, то предусматривается автоматический сдвигинтервала поиска корня вверх на величину равную b-a. В противном случаенаходится корень методом деления интервала пополам с вычислением функциив середине интервала. Процедура заканчивает работу, если значение функцииоказалось меньше произвольно заданного eps.

Величину следует выбиратьпримерно равной требуемой точности определения корня. При этом eps недолжно быть меньше, чем две единицы последнего разряда машинного слова,иначе возникнет зацикливание из-за округления при делении пополам.Хотя этот метод нулевого порядка и, следовательно, относится к самыммедленным методам, он применим к любой непрерывной функции. Тот факт,что эта процедура не содержит никакой проверки дифференцируемости, делаетее "надежной рабочей лошадкой" среди программ отыскания вещественныхкорней, которые предварительно были выделены.Произвольные переменные a и b являются (предположительно) концамиинтервала, внутри которого имеется нечетное число корней вещественнойфункции f[16,c.72].Function bisec(FexZ : TypeFun; a,b: realType;eps :RealType; var signal :boolean):RealType;Label 1,2,cont;var c, ya,yb,yc,st : RealType;beginsignal:=true;st:=b-a;ya:=FexZ(a);cont: yb:=FexZ(b);if ya*yb > 0 thenbeginya:=yb;a:=b;b:=b+st;goto cont;end;1: c:=(a+b)/2.0;yc:=FexZ(c);if (abs(a-b)< epsx) then goto 2;6if (yc*ya > 0) thenbegin a:=c; ya:=yc; goto 1 endelse begin b:=c; yb:=yc; goto 1 end;2: if (yc*ya > 0) thenbegin a:=c; ya:=yc endelse begin b:=c; yb:=yc end;bisec:=a+(a-b)*ya/(yb-ya);end;end;program testbisec;Uses Crt, BisecU, Kernel;varroot : RealType;signal : boolean;Function FexZer(x: realType): realType;beginFexZer:=sin(x)end;Beginroot:=bisec(FexZer,1.0, 4.0,0.1e-3,signal);write(root,' May be - ',pi);End.АЛГОРИТМ № 9б.

Интегрирование методом Рунге-КуттаПроцедура runge интегрирует системуyk' = f k ( x , y1, y2,..., yn )(k=1,2,..,n)обыкновенных дифференциальных уравнений методом Рунге-Кутта савтоматическим выбором шага.Параметры:x-начальное значение независимой переменной x;y[k]-начальные значения искомых функций y[k](x);n-порядок системы уравнений;f(x,y,n,z)-процедура, представляющая интегрированную систему[(совокупность функций f[k]), т.е. процедура, осуществляющая счет правых частейуравнений и выдающая в качестве результата массив z значений производныхy'[k] в любой заданной "точке" (x,y1,y2,..,[n]), т.е. z[k]=f[k](x,y[1],y[2],..,y[n]);eps,eta-значения, определяющие точность численного интегрирования;xfin-конец интервала интегрирования;yfin-выходной параметр, представляющий решение в точке x=xfin;prim-булевский параметр, значение которого при первом обращении впроцедуре runge должно быть true.

Если же нужно получить значения функций yв нескольких промежуточных точках x[0],x[1],..,x[n], то процедуру необходимовызвать повторно n раз с x=x[k], xfin=x[k+1] для k=0,1,..,n-1. В последнем случаеповторные обращения к runge для экономии машинного времени можноосуществлять с prim=false;comp-нелокализованный (глобальный) идентификатор вещественнойпроцедуры-функции comp (a,b,c), значением которой является абсолютное7значение разности мантисс a и b после того как порядки этих величинвыравнены до наибольшего порядка параметров a ,b и c ;rk1step(x,y,h,xh,yh)-процедура, интегрирующая систему уравнений на одномшаге h методом Рунге-Кутта с начальными значениями x,y[k].

Выходнымипараметрами являются xh=x+h и yh[k], причем последний вектор есть решениесистемы в точке xh. Параметры f,n,z входят в rk1step как нелокализованныеобъекты;ss,hs-глобальные параметры, значения которых не должны изменяться междуповторными обращениями к процедуре runge. При первом обращении к rungeзначения этих параметров могут быть не заданы. Переменная hs должна иметьтип real, a ss-тип integer.**Эта процедура, в отношении времени вычислений и ошибок округления,возможно, не будет оптимальной.{ Subroutine: Runge-Kutt 4 number from the "Communication of the ACM"algorithm 9 b , adapted by Kvashnin S.E. }function comp(a,b,c: RealType):RealType;varae,be,ce,i,aae : integer;co: RealType;function expon(x: RealType): integer;beginif x=0.0 then expon:=-37else expon:= round(0.4342944819*ln(abs(x))+1);end; { Expon }Beginae:=expon(a); be:=expon(b); ce:=expon(c);if ae < be then ae:=be;if ae < ce then ae:=ce;co:=1;aae:=abs(ae);if ae<>0 thenbeginfor i:=1 to aae do co:=co*10;if ae<0 then co:=1.0/co;end;comp:=abs(a-b)/co;End; { comp }UNIT RungeU;INTERFACEUSES Kernel;typeProcType = procedure ( T : RealType; var YI; n : integer; var YPI);var ss : integer;hs : RealType;8procedure Runge(Fex : ProcType; x: RealType; var y5; n :integer;var eps,eta,xfin : RealType; var prim : boolean;var yfin1);{Runge(FexRunge,x,y,n,eps,eta,xfin,prim,yfin); }Implementation{ Subroutine: Runge-Kutt 4 number from the "Communication of the ACM"algorithm 9 b , adapted by Kvashnin S.E.

}function comp(a,b,c: RealType) : RealType;var ae,be,ce,i,aae : longint;co,com: RealType;function expon(x: RealType): longint;beginif x=0.0 then expon:=-324else expon:= round(0.4342944819*ln(abs(x))+1);end; { Expon }Beginae:=expon(a); be:=expon(b); ce:=expon(c);if ae < be then ae:=be;if ae < ce then ae:=ce;co:=1;aae:=abs(ae);if ae<>0 thenbeginfor i:=1 to aae do co:=co*10;if ae<0 then co:=1.0/co;end;comp:=abs(a-b)/co;End; { comp }procedure Runge(Fex : ProcType; x: RealType; var y5; n :integer;var eps,eta,xfin : RealType; var prim : boolean;var yfin1);label 1,11,12,13,100;const nn=100; { max order of diff.equations }varx1,x2,x3,hk,joutz,y1,y2,y3yyfin: RealType;: integer;: boolean;: Array[1..nn] of realType;: RealTypeArray absolute y5;: RealTypeArray absolute yfin1;procedure rk1step(Fex : ProcType; x : RealType; var y1; h : RealType;9var xh : RealType; var yh1);varj,k: integer;a: array[1..5] of RealType;w: array[1..nn] of RealType;y : RealTypeArray absolute y1;yh: RealTypeArray absolute yh1;Begina[1]:=0.5*h;a[2]:=0.5*h;a[5]:=0.5*h;a[3]:=h; a[4]:=h; xh:=x;for k:=1 to n dobegin yh[k]:=y[k];w[k]:=y[k]; end;for j:=1 to 4 dobegin Fex(xh,w,n,z);xh:=x+a[j];for k:=1 to n dobegin yh[k]:=yh[k]+a[j+1]*z[k]/3;w[k]:=y[k]+a[j]*z[k];end; {k}end; {j}End; {rk1step}Begin1: if prim thenbeginh:=xfin-x; ss:=0;endelse h:=hs;out:=false;11: if((x+2.01*h-xfin)*h>0.0) thenbegin hs:=h; out:=true; h:=(xfin-x)/2.0; end;rk1step(Fex,x,y,2*h,x1,y1);12: rk1step(Fex,x,y,h,x2,y2);rk1step(Fex,x2,y2,h,x3,y3);for k:=1 to n doif comp(y1[k],y3[k],eta) > eps then goto 13;x:=x3;if out then goto 100;for k:=1 to n do y[k]:=y3[k];if ss=5 then begin ss:=0; h:=2*h end;ss:=ss+1; goto 11;13: h:=0.5*h; out:=false; x1:=x2;for k:=1 to n do y1[k]:=y2[k];goto 12;100:for k:=1 to n do yfin[k]:=y3[k];End; {Runge}{****************************************************************************}end.program testRunge;{---------------------------------------------------------------------------Subroutine RKF45 by H.A.WATTS, L.F.SHAMPINE , NEW Mexico10Adapted by S.Kvashnin, MHTS, January 1988----------------------------------------------------------------------------}program rungetest;uses Crt,rungeU,Kernel;const n=4;label 10,stop;type mass_n=array[1..n] of RealType;varecc,t,alfa,alfasq,eps,eta,tfinal,tprint,tout,h,savre,savae :realType;iflag,nfe,kop,init,jflag,kflag: integer;y,ypp,f1,f2,f3,f4,f5: mass_n;prim : boolean;procedure FexRunge(t: realType; var y1; n: integer; var yp1);varr: realType;w: RealTypeArray absolute Y1;yp: RealTypeArray absolute Yp1;begin{ test ORBIT }r:=w[1]*w[1]+w[2]*w[2];r:=r*sqrt(r)/alfasq;yp[1]:=w[3];yp[2]:=w[4];yp[3]:=-w[1]/r;yp[4]:=-w[2]/r;end;{---------------------------------------------------------------------------}beginalfa:=pi*0.25;alfasq:=alfa*alfa;ecc:=0.25;t:=0;y[1]:=1.0-ecc;y[2]:=0;y[3]:=0;y[4]:=alfa*sqrt((1.0+ecc)/(1.0-ecc));eps:=1e-9;eta:=0.0;tfinal:=12.0;tprint:=0.5;prim:=true;tout:=t;10: Runge(FexRunge,t,y,n,eps,eta,tout, prim,ypp);writeln(iflag,t,y[1],y[2]);goto 10;stop:end.