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

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

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

右表は塩素 Cl₂ の水への溶解度のデータである．yは Cl₂の気相のモル分率，xはこれに平衡な水中の濃度（モル分率）である．Cl₂は水中で解離するので，ヘンリー法則には従わない．

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

のように解離して，水溶液中には塩素が分子Cl₂，イオンCl^-，HOCl　の３つの形で存在する．水中の塩素全濃度xは，

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

なので，

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

この２式より，
x=(1/H)y+((K/H)y)^1/3
となる．すなわち，Cl₂の溶解度データは，
x=ay+by^1/3
の関係と予想される．

実習の手順

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

相関式： x=1.23x10^-3y+5.05x10^-4y^1/3

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

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

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

▲目次