1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
|
PROCEDURE mtt_implicit(VAR xnew,x,dx : StateVector;
VAR AA : StateMatrix;
VAR AAx : StateVector;
dt : REAL;
Nx : INTEGER;
VAR open : StateVector);
VAR
i,ii,j,jj : INTEGER;
BB,xsub : StateVector;
AAsub : StateMatrix;
BEGIN{mtt_implicit}
ii := 0;
FOR i := 1 TO Nx DO {BB is (1-A*dt)*x +dx*dt}
BEGIN
IF open[i]<0.5 THEN
BEGIN
ii := ii+1; jj := 0;
BB[ii] := AAx[i] + dt*dx[i];
FOR j := 1 TO Nx DO
IF open[j]<0.5 THEN
BEGIN
jj := jj+1;
AAsub[ii,jj] := AA[i,j];
END;
END;
END;
Nx := ii; {New size}
mtt_solve(xsub,AAsub,BB,Nx); {Solve AAx=BB}
ii := 0;
FOR i := 1 TO Nx DO
IF open[i]<0.5 THEN
BEGIN
ii := ii+1;
xnew[i] := xsub[ii];
END ELSE
|
|
>
|
|
|
|
|
|
|
|
>
|
|
|
|
>
|
>
>
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
PROCEDURE mtt_implicit(VAR xnew : StateVector;
x,dx : StateVector;
AA : StateMatrix;
AAx : StateVector;
dt : REAL;
Nx : INTEGER;
open : StateVector);
VAR
i,ii,j,jj,Nxsub : INTEGER;
BB,xsub : StateVector;
AAsub : StateMatrix;
BEGIN{mtt_implicit}
ii := 0;
FOR i := 1 TO Nx DO {BB is (1-A*dt)*x +dx*dt}
BEGIN
IF open[i]<0.5 THEN
BEGIN
ii := ii+1; jj := 0;
BB[ii] := AAx[i] + dt*dx[i];
FOR j := 1 TO Nx DO
IF open[j]<0.5 THEN
BEGIN
jj := jj+1;
AAsub[ii,jj] := AA[i,j];
END;
END;
END;
Nxsub := ii; {New size}
mtt_solve(xsub,AAsub,BB,Nxsub); {Solve AAx=BB}
ii := 0;
FOR i := 1 TO Nx DO
BEGIN
IF open[i]<0.5 THEN
BEGIN
ii := ii+1;
xnew[i] := xsub[ii];
END
ELSE
xnew[i] := 0.0;
writeln(i,ii);
END;
END;{mtt_implicit}
|