Var I : integer;

sum_x,sum_y,sum_xy,sum_x2,

sum_y2,xi,yi,sxy,sxx,

syy : real;

begin { linfit }

sum_x:=0.0;

sum_y:=0.0;

sum_xy:=0.0;

sum_x2:=0.0;

sum_y2:=0.0;

for i:=1 to n do

begin

xi:=x[i];

yi:=y[i];

sum_x:=sum_x+xi;

sum_y:=sum_y+yi;

sum_xy:=sum_xy+xi*yi;

sum_x2:=sum_x2+xi*xi;

sum_y2:=sum_y2+yi*yi;

end;

sxx:=sum_x2-sum_x*sum_x/n;

sxy:=sum_xy-sum_x*sum_y/n;

syy:=sum_y2-sum_y*sum_y/n;

b:=sxy/sxx;

a:=((sum_x2*sum_y-sum_x*sum_xy)/n)/sxx;

for i:=1 to n do

y_calc[i]:=a+b*x[i]

end; { LINFIT }

Программа 16. Приближенная линеаризация опытных данных (вар.2).

procedure linfit2(x,y: ary;

var y_calc: ary;

var a,b: real;

n: integer);

{ fit a straight line (y_calc) through n sets of x and y pairs of points }

Var I : integer;

sum_x,sum_y,sum_xy,sum_x2,

sum_y2,xi,yi,sxy,syy,

sxx : real;

begin { linfit }

sum_x:=0.0;

sum_y:=0.0;

sum_xy:=0.0;

sum_x2:=0.0;

sum_y2:=0.0;

for i:=1 to n do

begin

xi:=x[i];

yi:=y[i];

sum_x:=sum_x+xi;

sum_y:=sum_y+yi;

sum_xy:=sum_xy+xi*yi;

sum_x2:=sum_x2+xi*xi;

sum_y2:=sum_y2+yi*yi;

end;

sxx:=sum_x2-sum_x*sum_x/n;

sxy:=sum_xy-sum_x*sum_y/n;

syy:=sum_y2-sum_y*sum_y/n;

b:=sxy/sxx;

a:=((sum_x2*sum_y-sum_x*sum_xy)/n)/sxx;

correl_coef:=sxy/sqrt(sxx*syy);

see:=sqrt((sum_y2-a*sum_y-b*sum_xy)/(n-2));

sigma_b:=see/sqrt(sxx);

sigma_a:=sigma_b*sqrt(sum_x2/n);

for i:=1 to n do

y_calc[i]:=a+b*x[i]

end; { LINFIT }

Программа 17. Задача линейного приближения с помощью МНК (вар.2).

program cfit1A;

{ Pascal program to perform a linear least-squares fit

with random generator }

const max = 20;

type index = 1..max;

ary = array[index] of real;

var x,y,y_calc : ary;

n : integer;

first,done : boolean;

seed,a,b : real;

function random(dummy: integer): real;

{ random number 0-1 }

{ define seed=4.0 as global }

const pi = 3.14159;

var x : real;

I : integer;

begin { RANDOM }

x:=seed+pi;

x:=exp(5.0*ln(x));

seed:=x-trunc(x);

random:=seed

end; { RANDOM }

procedure get_data(var x,y: ary;

var n: integer);

{ get values for n and arrays x,y }

{ y is randomly scattered about a straight line }

const a = 2.0;

b = 5.0;

var i,j : integer;

fudge : real;

begin

write('Fudge? (<0 to terminate) ');

readln(fudge);

if fudge<0.0 then done:=true

else

begin

repeat

write('How many points? ');

readln(n)

until (n>2) and (n<=max);

if first then first:=false else ClrScr;

for i:=1 to n do

begin

j:=n+1-i;

x[i]:=j;

y[i]:=(a+b*j)*(1.0+(2.0*random(0)-1.0)*fudge)

end { for-loop }

end { if }

end; { procedure get_data }

procedure write_data;

{ print out the answers }

Var I : integer;

begin

writeln;

writeln(' I X Y');

for i:=1 to n do

writeln(i:3,x[i]:8:1,y[i]:9:2);

writeln

end; { write_data }

begin { MAIN program }

ClrScr;

seed:=4.0;

first:=true;

done:=false;

repeat

get_data(x,y,n);

if not done then

begin

write_data;

{ ***** ---> more lines to be added here ********* }

end

until done

end.

Программа 18. Решение системы уравнений методом Гаусса-Жордана.

(Требуется согласовать упрощения)

procedure gaussj

(var b: ary2s; { square matrix of coefficients }

y: arys; { constant vector }

var coef: arys; { solution vector }

ncol: integer; { order of matrix }

var error: boolean); { true if matrix singular }

{ Gauss Jordan matrix inversion and solution }

{ B(n,n) coefficient matrix becomes inverse }

{ Y(n) original constant vector }

{ W(n,m) constant vector(s) become solution vector }

{ DETERM is the determinant }

{ ERROR=1 if singular }

{ INDEX(n,3) }

{ NV is number of constant vectors }

label 99;

var

w : array[1..maxc,1..maxc] of real;

index : array[1..maxc,1..3] of integer;

i,j,k,l,nv,

irow,icol,

n,l1 : integer;

determ,pivot,

hold,sum,t,

ab,big : real;

procedure swap(var a,b: real);

var hold : real;

begin { swap }

hold:=a;

a:=b;

b:=hold

end; { procedure swap }

procedure gausj2;

label 98;

var i,j,k,l,l1 : integer;

procedure gausj3;

var l : integer;

begin { procedure gausj3 }

{ interchange rows to put pivot on diagonal }

if irow<>icol then

begin

determ:=-determ;

for l:=1 to n do

swap(b[irow,l],b[icol,l]);

if nv>0 then

for l:=1 to nv do

swap(w[irow,l],w[icol,l])

end { if iroe<>icol }

end; { gausj3 }

begin { procedure gausj2 }

{ actual start of gaussj }

error:=false;

nv:=1; { single constant vector }

n:=ncol;

for i:=1 to n do

begin

w[i,1]:=y[i]; { copy constant vector }

index[i,3]:=0

end;

determ:=1.0;

for i:=1 to n do

begin

{ search for largest element }

big:=0.0;

for j:=1 to n do

begin

if index[j,3]<>1 then

begin

for k:=1 to n do

begin

if index[k,3]>1 then

begin

writeln('ERROR: matrix is singular');

error:=true;

goto 98 { abort }

end;

if index[k,3]<1 then

if abs(b[j,k])>big then

begin

irow:=j;

icol:=k;

big:=abs(b[j,k])

end

end { k-loop }

end

end; { j-loop }

index[icol,3]:=index[icol,3]+1;

index[i,1]:=irow;

index[i,2]:=icol;

gausj3; { further subdivision of gaussj }

{ divide pivot row by pivot column }

pivot:=b[icol,icol];

determ:=determ*pivot;

b[icol,icol]:=1.0;

for l:=1 to n do

b[icol,l]:=b[icol,l]/pivot;

if nv>0 then

for l:=1 to nv do

w[icol,l]:=w[icol,l]/pivot;

{ reduce nonpivot rows }

for l1:=1 to n do

begin

if l1<>icol then

begin

t:=b[l1,icol];

b[l1,icol]:=0.0;

for l:=1 to n do

b[l1,l]:=b[l1,l]-b[icol,l]*t;

if nv>0 then

for l:=1 to nv do

w[l1,l]:=w[l1,l]-w[icol,l]*t;

end { if l1<>icol }

end

end; { i-loop }

98:

end; { gausj2 }

begin { gaus-jordan main program }

gausj2; { first half of gaussj }

if error then goto 99;

{ interchange columns }

for i:=1 to n do

begin

l:=n-i+1;

if index[l,1]<>index[l,2] then

begin

irow:=index[l,1];

icol:=index[l,2];

for k:=1 to n do

swap(b[k,irow],b[k,icol])

end { if index }

end; { i-loop }

for k:=1 to n do

if index[k,3]<>1 then

begin

writeln(chr(7),'ERROR: matrix is singular');

error:=true;

goto 99 { abort }

end;

for i:=1 to n do

coef[i]:=w[i,1];

99:

end; { procedure gaussj }

Программа 19. Вычисление функции ошибок распределения Гаусса (вар.1).

program erfd4;

{ evaluation of the gaussian error function }

var x,er,ec : real;

done : boolean;

function erf(x: real): real;

{ infinite series expansion of the Gaussian error function }

const sqrtpi = 1.7724538;

t2 = 0.66666667;

t3 = 0.66666667;

t4 = 0.07619048;

t5 = 0.01693122;

t6 = 3.078403E-3;

t7 = 4.736005E-4;

t8 = 6.314673E-5;

t9 = 7.429027E-6;

t10 = 7.820028E-7;

t11 = 7.447646E-8;

t12 = 6.476214E-9;

var x2,sum : real;