?? unit2.pas
字號:
unit Unit2;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;
Function GASDEV:real;
Procedure SVDVAR(V:matrx2; MA, NP:integer; W:array of real;
var CVM:matrx2; NCVM:integer);
Procedure SVDFIT(X, Y, SIG:array of real; NDATA:integer;
var A:array of real; MA:integer;var U, V:matrx2;var W:array of real;
MP, NP:integer;var CHISQ:real; FUNCS:string);
implementation
Function GASDEV:real;
var
V1,V2,FAC,R:real;
begin
If ISET^= 0 Then
begin
repeat
V1:=2 * Random - 1;
V2:=2 * Random - 1;
R:=Sqr(V1) + Sqr(V2);
until (R < 1);
FAC:=Sqrt(-2 * Ln(R) / R);
GSET^:=V1 * FAC;
GASDEV:=V2 * FAC;
ISET^:=1;
end
Else
begin
GASDEV:=GSET^;
ISET^:=0;
end;
end;
procedure SVDCMP(A:matrx2; M, N:integer;var W:array of real; V:matrx2);
label 1,2,3;
var
RV1:array [0..100] of real;
I,J,L,ITS,JJ,NM,K:integer;
G,S,SCALE1,C,F,H,X,Y,Z,Sgn,ANORM,AAAAA:real;
begin
If M < N Then
ShowMessage('You must augment A with extra zero rows.');
G:=0;
SCALE1:=0;
ANORM:=0;
For I:=1 To N do
begin
L:=I + 1;
RV1[I]:=SCALE1 * G;
G:=0;
S:=0;
SCALE1:=0;
If I <= M Then
begin
For K:=I To M do
SCALE1:=SCALE1 + Abs(A[K, I]);
If SCALE1 <> 0 Then
begin
For K:=I To M do
begin
A[K, I]:=A[K, I] / SCALE1;
S:=S + A[K, I] * A[K, I];
end;
F:=A[I, I];
If F > 0 Then
SGN:=1
else
SGN:=-1;
G:=-Sqrt(S) * Sgn;
H:=F * G - S;
A[I, I]:=F - G;
If I <> N Then
begin
For J:=L To N DO
begin
S:=0;
For K:=I To M do
S:=S + A[K, I] * A[K, J];
F:=S / H;
For K:=I To M do
A[K, J]:=A[K, J] + F * A[K, I];
end;
end;
For K:=I To M do
A[K, I]:=SCALE1 * A[K, I];
end;
end;
W[I]:=SCALE1 * G;
G:=0;
S:=0;
SCALE1:=0;
If (I <= M) And (I <> N) Then
begin
For K:=L To N do
SCALE1:=SCALE1 + Abs(A[I, K]);
If SCALE1 <> 0 Then
begin
For K:=L To N do
begin
A[I, K]:=A[I, K] / SCALE1;
S:=S + A[I, K] * A[I, K];
end;
F:=A[I, L];
If F > 0 Then
SGN:=1
else
SGN:=-1;
G:=-Sqrt(S) * Sgn;
H:=F * G - S;
A[I, L]:=F - G;
For K:=L To N do
RV1[K]:=A[I, K] / H;
If I <> M Then
begin
For J:=L To M do
begin
S:=0;
For K:=L To N do
S:=S + A[J, K] * A[I, K];
For K:=L To N do
A[J, K]:=A[J, K] + S * RV1[K];
end;
end;
For K:=L To N do
A[I, K]:=SCALE1 * A[I, K];
end;
end;
If ANORM > Abs(W[I]) + Abs(RV1[I]) Then
ANORM:=ANORM
Else
ANORM:=Abs(W[I]) + Abs(RV1[I]);
end;
For I:=N DownTo 1 do
begin
If I < N Then
begin
If G <> 0 Then
begin
For J:=L To N do
V[J, I]:=(A[I, J] / A[I, L]) / G;
For J:=L To N do
begin
S:=0;
For K:=L To N do
S:=S + A[I, K] * V[K, J];
For K:=L To N do
V[K, J]:=V[K, J] + S * V[K, I];
end;
end;
For J:=L To N do
begin
V[I, J]:=0;
V[J, I]:=0;
end;
end;
V[I, I]:=1;
G:=RV1[I];
L:=I;
end;
For I:=N DownTo 1 do
begin
L:=I + 1;
G:=W[I];
If I < N Then
begin
For J:=L To N do
A[I, J]:=0;
end;
If G <> 0 Then
begin
G:=1 / G;
If I <> N Then
begin
For J:=L To N do
begin
S:=0;
For K:=L To M do
S:=S + A[K, I] * A[K, J];
F:=(S / A[I, I]) * G;
For K:=I To M do
A[K, J]:=A[K, J] + F * A[K, I];
end;
end;
For J:=I To M do
begin
A[J, I]:=A[J, I] * G;
end;
end
Else
For J:=I To M do
A[J, I]:=0;
A[I, I]:=A[I, I] + 1;
end;
For K:=N DownTo 1 do
begin
For ITS:=1 To 30 do
begin
For L:=K DownTo 1 do
begin
NM:=L - 1;
If Abs(RV1[L]) + ANORM = ANORM Then GoTo 2;
If Abs(W[NM]) + ANORM = ANORM Then GoTo 1;
end;
1: C:=0;
S:=1;
For I:=L To K do
begin
F:=S * RV1[I];
If Abs(F) + ANORM <> ANORM Then
begin
G:=W[I];
H:=Sqrt(F * F + G * G);
W[I]:=H;
H:=1 / H;
C:=(G * H);
S:=-(F * H);
For J:=1 To M do
begin
Y:=A[J, NM];
Z:=A[J, I];
A[J, NM]:=(Y * C) + (Z * S);
A[J, I]:=-(Y * S) + (Z * C);
end;
end;
end;
2: Z:=W[K];
If L = K Then
begin
If Z < 0 Then
begin
W[K]:=-Z;
For J:=1 To N do
V[J, K]:=-V[J, K];
end;
GoTo 3;
end;
If ITS = 30 Then
ShowMessage('No convergence in 30 iterations');
X:=W[L];
NM:=K - 1;
Y:=W[NM];
G:=RV1[NM];
H:=RV1[K];
F:=((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2 * H * Y);
G:=Sqrt(F * F + 1);
If F > 0 Then
Sgn:=1
else
Sgn:=-1;
F:=((X - Z) * (X + Z) + H * ((Y / (F + Abs(G) * Sgn)) - H)) / X;
C:=1;
S:=1;
For J:=L To NM do
begin
I:=J + 1;
G:=RV1[I];
Y:=W[I];
H:=S * G;
G:=G * C;
Z:=Sqrt(F * F + H * H);
RV1[J]:=Z;
C:=F / Z;
S:=H / Z;
F:=(X * C) + (G * S);
G:=-(X * S) + (G * C);
H:=Y * S;
Y:=Y * C;
For NM:=1 To N do
begin
X:=V[NM, J];
Z:=V[NM, I];
V[NM, J]:=(X * C) + (Z * S);
V[NM, I]:=-(X * S) + (Z * C);
end;
Z:=Sqrt(F * F + H * H);
W[J]:=Z;
If Z <> 0 Then
begin
Z:=1 / Z;
C:=F * Z;
S:=H * Z;
end;
F:=(C * G) + (S * Y);
X:=-(S * G) + (C * Y);
For NM:=1 To M do
begin
Y:=A[NM, J];
Z:=A[NM, I];
A[NM, J]:=(Y * C) + (Z * S);
A[NM, I]:=-(Y * S) + (Z * C);
end;
end;
RV1[L]:=0;
RV1[K]:=F;
W[K]:=X;
end;
3: AAAAA:=1;
end;
end;
procedure SVBKSB(U:matrx2; W:array of real; V:matrx2;
M, N:integer; B:array of real;var X:array of reaL);
var
TMP:array[0..100] of real;
I,J,JJ:integer;
S:real;
begin
For J:=1 To N do
begin
S:=0;
If W[J] <> 0 Then
begin
For I:=1 To M do
S:=S + U[I, J] * B[I];
S:=S / W[J];
end;
TMP[J]:=S;
end;
For J:=1 To N do
begin
S:=0;
For JJ:=1 To N do
S:=S + V[J, JJ] * TMP[JJ];
X[J]:=S;
end;
end;
Procedure FPOLY(X:real; var P:array of real; NP:integer);
var
J:integer;
begin
P[1]:=1;
For J:=2 To NP do
P[J]:=P[J - 1] * X;
end;
Procedure FLEG(X:real; var PL:array of real; NL:integer);
var
J:integer; TWOX,F1,F2,D:real;
begin
PL[1]:=1;
PL[2]:=X;
If NL > 2 Then
begin
TWOX:=2 * X;
F2:=X;
D:=1;
For J:=3 To NL do
begin
F1:=D;
F2:=F2 + TWOX;
D:=D + 1;
PL[J]:=(F2 * PL[J - 1] - F1 * PL[J - 2]) / D;
end;
end;
end;
Procedure SVDVAR(V:matrx2; MA, NP:integer; W:array of real;
var CVM:matrx2; NCVM:integer);
var
WTI:array[0..20] of real;
I,K,J:integer; SUM1:real;
begin
For I:=1 To MA do
begin
WTI[I]:=0 ;
If W[I] <> 0 Then WTI[I]:=1 / (W[I] * W[I]);
end;
For I:=1 To MA do
begin
For J:=1 To I do
begin
SUM1:=0;
For K:=1 To MA do
SUM1:=SUM1 + V[I, K] * V[J, K] * WTI[K];
CVM[I, J]:=SUM1;
CVM[J, I]:=SUM1;
end;
end;
end;
Procedure SVDFIT(X, Y, SIG:array of real; NDATA:integer;
var A:array of real; MA:integer;var U, V:matrx2;var W:array of real;
MP, NP:integer;var CHISQ:real; FUNCS:string);
var
B:array[0..1000] of real; AFUNC:array[0..50] of real;
I,J:integer; TMP,WMAX,THRESH,SUM1:real;
const
TOL=0.00001;
begin
For I:=1 To NDATA do
begin
If FUNCS = 'FPOLY' Then FPOLY(X[I], AFUNC, MA);
If FUNCS = 'FLEG' Then FLEG(X[I], AFUNC, MA);
TMP:=1 / SIG[I];
For J:=1 To MA do
U[I, J]:=AFUNC[J] * TMP;
B[I]:=Y[I] * TMP;
end;
SVDCMP(U, NDATA, MA, W, V);
WMAX:=0;
For J:=1 To MA do
If W[J] > WMAX Then WMAX:=W[J];
THRESH:=TOL * WMAX;
For J:=1 To MA do
If W[J] < THRESH Then W[J]:=0;
SVBKSB(U, W, V, NDATA, MA, B, A);
CHISQ:=0;
For I:=1 To NDATA do
begin
If FUNCS = 'FPOLY' Then FPOLY(X[I], AFUNC, MA);
If FUNCS = 'FLEG' Then FLEG(X[I], AFUNC, MA);
SUM1:=0;
For J:= 1 To MA do
SUM1:=SUM1 + A[J] * AFUNC[J];
CHISQ:=CHISQ + Sqr((Y[I] - SUM1) / SIG[I]);
end;
end;
end.
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -