平衡関係の相関 −パラメータ推定−

 物性値や装置の特性などのデータを解析する場合には必ず何らかの理論が想定され,その理論により導かれた理論式の形にもとずいてデータが相関される.化学工学における理論式はほとんど非線形なので,データの相関には一般の1次式,2次式などの最小2乗法は適用できない.
 ここでは,データ→これを相関する理論式→非線形最小2乗法による相関→データを表わす相関式の作成,の手順を実習する.

演習例題) 塩素の水への溶解平衡関係

右表は塩素 Cl2 の水への溶解度のデータである.yは Cl2の気相のモル分率,xはこれに平衡な水中の濃度(モル分率)である.Cl2は水中で解離するので,ヘンリー法則には従わない.

x
(20℃)
x
(40℃)
y
0.000111 0.0001056 0.00658
0.000146 0.0001362 0.0131
0.000237 0.000210 0.0394
0.000308 0.000263 0.0657
0.000452 0.000363 0.131
0.000578 0.000448 0.197
0.000696 0.000524 0.263
0.000812 0.000599 0.328
0.000923 0.000668 0.394
0.001033 0.000732 0.460

 このことは次の考察で説明される.水に溶解したCl2の一部は,

のように解離して,水溶液中には塩素が分子Cl2,イオンCl-,HOCl の3つの形で存在する.水中の塩素全濃度xは,

である.なお,である.この反応の平衡定数をKとすると,

なので,

である.分子のままで水中に存在する塩素の濃度 xCl2についてはヘンリーの法則が成り立つ.
y=HxCl2

この2式より,
x=(1/H)y+((K/H)y)1/3
となる.すなわち,Cl2の溶解度データは,
x=ay+by1/3
の関係と予想される.

実習の手順

@ 方眼紙に横軸 x,縦軸 y の座標上に表の40℃,20℃のデータをプロットする.
A この値より,学籍番号末尾が0と5,1と6,2と7,3と8,4と9,毎にそれぞれ15℃,18℃,25℃,28℃,32℃の値を6点推定してデータとしてプロットする.
B この6つのデータを数値化して表とする.
C 非線形最小2乗法プログラムにより,このデータを式:x =ay+by1/3 で相関して,パラメータ a,b を求める.
(x,yの関係が通常と逆で,yのほうから計算することに注意)
D 相関式をグラフ上に描く.

相関式: x=1.23x10-3y+5.05x10-4y1/3

非線形最小2乗法プログラム

(プログラムの解説は省略する.基本的には第2回の非線形連立方程式の解法と同じである.ただし,式の数が未知数の数より多い.)

10 ' ヒセンケイサイショウニジョウホウ (トネ カオル "プレイマイコンシリーズ1 BASIC" バイフウカン 1981)
12 READ N, M
13 IMAX = 50 * N * M: N1 = N + 1: N2 = N + 2: M1 = M + N + 1: M2 = M1 + 1: M3 = M2 + 1
16 DIM x(N), F(M), B(M3, N2)
102 '
103 ' FUNCTION (Cl2/Water キエキヘイコウ)
107 G$ = "XX=X(1)*y+X(2)*y^(1/3)"
110 GOTO 8000
120 '
125 ' FUNCTION VALUES
150 FOR I = 1 TO M: READ Y, XX
151 F(I) = XX - (x(1) * Y + x(2) * Y ^ (1 / 3))
152 NEXT I: RESTORE 170: RETURN
155 ' X(1),X(2)..ノ カズ, データ クミスウ , X(1),X(2)..ノ ショキチ
160 DATA 2, 6, 0.0001, 0.003
165 'データ y vs. xx
170 DATA 0.01, 0.0001
171 DATA 0.05,0.00025
172 DATA 0.14,0.00044
173 DATA 0.23,0.0006
174 DATA 0.34,0.00078
175 DATA 0.43,0.0009
1000 '-------JACOBIAN
1010 FOR I = 1 TO M: B(I, N1) = F(I): NEXT I
1020 IF JS = 1 THEN 1070: 'CENTRAL DIFF
1030 FOR I = 1 TO M: B(I, N2) = F(I): NEXT I
1040 FOR J = 1 TO N: x(J) = x(J) + H: GOSUB 120: NF = NF + 1
1050 FOR I = 1 TO M: B(I, J) = (F(I) - B(I, N1)) / H: NEXT I
1060 x(J) = x(J) - H: NEXT J: RETURN
1070 FOR J = 1 TO N: x(J) = x(J) - H: GOSUB 120: NF = NF + 1
1080 FOR I = 1 TO M: B(I, N + 2) = F(I): NEXT I
1090 x(J) = x(J) + 2 * H: GOSUB 120: NF = NF + 1
1100 FOR I = 1 TO M: B(I, J) = (F(I) - B(I, N2)) / (2 * H): NEXT I
1110 x(J) = x(J) - H: NEXT J
1120 FOR I = 1 TO M: F(I) = B(I, N1): NEXT I
1130 RETURN
2000 '-------DELTA-X
2010 B(N, N2) = B(N, N1) / B(M1, N): IF N = 1 THEN 2060
2020 FOR K = N - 1 TO 1 STEP -1: B(K, N2) = B(K, N1)
2030 FOR J = K + 1 TO N: B(K, N2) = B(K, N2) - B(K, J) * B(J, N2): NEXT J
2040 B(K, N2) = B(K, N2) / B(M1, K)
2050 NEXT K
2060 FOR I = 1 TO N: B(M3, I) = x(I) - B(I, N2): NEXT I
2070 FOR I = 1 TO N: B(M2, I) = x(I): x(I) = B(M3, I): NEXT I
2080 GOSUB 120: NF = NF + 1
2090 FOR I = 1 TO N: x(I) = B(M2, I): NEXT I
2100 RETURN
3000 '----TEST DELTA-X
3010 SN = 0: FOR I = 1 TO M: SN = SN + F(I) ^ 2: NEXT I
3020 IF SUMSQ < SN THEN SW = 1: RETURN
3030 SW = 0: RETURN
4000 '-------REFACTOROZATION
4010 RR = RINCR * RINCR - 1: V2 = V * V: IF IC > 1 THEN 4070
4020 FOR I = 1 TO N: B(I, I) = B(M1, I): NEXT I
4030 FOR I = 1 TO N: FOR J = N1 TO I STEP -1
4040 S = 0: FOR K = 1 TO I: S = S + B(K, I) * B(K, J): NEXT K
4050 B(J + 1, I) = S
4060 NEXT J: NEXT I
4070 RV = RR * V2: FOR I = 1 TO N: B(I + 1, I) = B(I + 1, I) + RV: NEXT I
4080 FOR K = 1 TO N: S = B(K + 1, K): IF K = 1 THEN 4100
4090 FOR I = 1 TO K - 1: S = S - B(I, K) * B(I, K): NEXT I
4100 B(K, K) = SQR(S): B(M1, K) = B(K, K) * SGN(B(M1, K))
4110 FOR J = K + 1 TO N1: S = B(J + 1, K): IF K = 1 THEN 4130
4120 FOR I = 1 TO K - 1: S = S - B(I, K) * B(I, J): NEXT I
4130 B(K, J) = S / B(M1, K)
4140 NEXT J: NEXT K
4150 RETURN
5000 '---------HOUSEHOLDER
5010 FOR K = 1 TO N: S = 0
5020 FOR I = K TO M + K: S = S + B(I, K) ^ 2: NEXT I
5030 S = SQR(S): BT = 1 / (S * (S + ABS(B(K, K))))
5040 IF B(K, K) = 0 THEN B(M1, K) = -S: B(K, K) = S: GOTO 5060
5050 B(M1, K) = -SGN(B(K, K)) * S: B(K, K) = B(K, K) + SGN(B(K, K)) * S
5060 FOR J = K + 1 TO N1: S = 0
5070 FOR L = K TO M + K: S = S + B(L, K) * B(L, J): NEXT L
5080 S = S * BT
5090 FOR I = K TO M + K: B(I, J) = B(I, J) - S * B(I, K): NEXT I
5100 NEXT J: NEXT K
5110 RETURN
8000 '---------- NONLINEAR LEAST SQUARES
8020 PRINT : PRINT "---------"
8030 PRINT "INITIAL CONDITIONS": PRINT
8040 PRINT "NO. OF VARIABLES= "; N: PRINT "NO. OF FUNCTIONS= "; M
8050 IMAX = 50 * N * M: N1 = N + 1: N2 = N + 2: M1 = M + N + 1: M2 = M1 + 1: M3 = M2 + 1
8060 PRINT "MAX NO. OF FUNCTION EVALUATIONS= "; IMAX
8070 '
8090 EPS = .00001: H = .00001
8100 PRINT "STEP SIZE IN NUMERICAL DIFFERENTIATION": PRINT "H="; H
8110 PRINT "CONVERGENCE CRITERION": PRINT "EPSILON="; EPS: PRINT
8120 NIT = 0: RINCR = 1.5: RDECR = .5: NF = 0: ' NO. OF F-EVALUATIONS
8130 FOR I = 1 TO N: READ x(I): PRINT "x("; I; ")="; x(I): NEXT I
8140 '
8143 PRINT G$
8146 PRINT : PRINT "INITIAL VALUE"
8150 FOR I = 1 TO N: PRINT USING "X(#)="; I; : PRINT x(I); : NEXT I: PRINT
8160 NIT = NIT + 1
8180 GOSUB 120: ' F(I,X )
8190 NF = NF + 1: S = 0
8200 FOR I = 1 TO M: S = S + F(I) ^ 2: NEXT I: SUMSQ = S
8210 PRINT : FOR I = 1 TO M: PRINT "F("; I; ")="; F(I): B(I, N1) = F(I): NEXT I
8220 PRINT : PRINT "SUMSQ="; SUMSQ
8230 GOSUB 1000: ' JACOBIAN
8240 IF NIT > 1 THEN 8270
8250 S = 0: FOR I = 1 TO M: FOR J = 1 TO N: S = S + B(I, J) ^ 2: NEXT J, I
8260 V = SQR(S / (N * M))
8270 IC = O: ' INNER COUNTER
8280 FOR I = 1 TO N: FOR J = 1 TO N1: B(M + I, J) = 0: NEXT J
8285 B(M + I, I) = V: NEXT I
8290 GOSUB 5000: ' HOUSEHOLDER
8300 GOSUB 2000: ' DX
8310 GOSUB 3000: ' TEST DX
8320 '------------TEST CONVERGENCE
8330 P = ABS(x(1)): PP = ABS(B(1, N2)): IF N = 1 THEN 8360
8340 FOR I = 2 TO N: IF ABS(x(I)) > P THEN P = ABS(x(I))
8350 IF ABS(B(I, N2)) > PP THEN PP = ABS(B(I, N2))
8355 NEXT I
8360 IF P < 1 THEN P = 1
8363 IF PP < EPS * P THEN 8510
8366 IF SW = 0 THEN 8390
8370 IC = IC + 1: GOSUB 4000
8380 V = RINCR * V: IF V > 1000000! THEN 8610
8385 GOTO 8300
8390 JS = 0: IF PP < EPS * P * 100 THEN JS = 1
8400 FOR I = 1 TO N: x(I) = B(M3, I): NEXT I
8410 SUMSQ = SN
8430 PRINT "NO. OF ITERATIONS="; NIT;
8440 FOR I = 1 TO N: PRINT "X("; I; " )="; x(I); : NEXT I
8450 IF IC > 0 THEN 8460
8455 V = V * RDECR
8460 IF NF > IMAX THEN 8620
8465 NIT = NIT + 1
8470 FOR I = 1 TO M: B(I, N1) = F(I): NEXT I
8480 GOTO 8220
8490 '--------- SOLUTION
8510 PRINT "サイテキチ";
8520 IF SW = 1 THEN SN = SUMSQ: GOTO 8540
8530 FOR I = 1 TO N: x(I) = B(M3, I): NEXT I
8540 FOR I = 1 TO N: PRINT " X("; I; ")="; x(I); : NEXT I
8548 RESTORE 170
8549 PRINT "No.,y(データ),x(データ),x(アテハメチ)";
8550 FOR I = 1 TO M: READ Y, XX
8551 PRINT "No."; I; Y; XX; XX - F(I);
8552 NEXT I
8560 'PRINT "SUN OF SQUARES="; SN: PRINT "NO. OF ITERATIONS="; NIT
8570 'PRINT "NO. OF FUNCTION EVALUATIONS="; NF
8590 END
8600 PRINT "TOO SMALL EPSR": STOP
8610 PRINT "TOO BIG V": STOP
8620 PRINT "NO CONVERGENCE WITHIN "; NF; " FUNCTION CALLS": STOP


▲目次