を解くことは，f(x)→0となるx=(x1,x2,････xn)T を求めることである。
一変数の非線形方程式 f(x)=0 の解法であるNewton-Raphson法：
xk+1=xk- f(xk)/f'(xk) (k=0,1,2････)
方程式を解く－非線形方程式－
連立非線形方程式の数値解法
n変数の連立非線形方程式
を解くことは，f(x)→0となるx=(x1,x2,････xn)T を求めることである。
一変数の非線形方程式
f(x)=0 の解法であるNewton-Raphson法：
xk+1=xk-
f(xk)/f'(xk) (k=0,1,2････)
が
n 元連立非線形方程式に対し拡張できる。この疑似Newton法による非線形連立方程式解法の計算手順は以下のようである。
1 (ﾌﾟﾛｸﾞﾗﾑ350-370) 適当な初期ﾍﾞｸﾄﾙ
x(0)を与える。k=0,1,2････に対し以下を反復計算する。
2 (ﾌﾟﾛｸﾞﾗﾑ380-390)
f(x(k))の値を計算する。各ﾍﾞｸﾄﾙの値は0でない。
3 (ﾌﾟﾛｸﾞﾗﾑ30-440,ｻﾌﾞﾙｰﾁﾝ3000) ∂fi
/∂xj を前進差分で近似して求める。
これをJacobi(ﾔｺﾋﾞ)行列と呼ぶ。
4 (ﾌﾟﾛｸﾞﾗﾑ450-550,ｻﾌﾞﾙｰﾁﾝ1000,2000) 次の連立１次方程式をGaussの消去法で解き，解Δx(k)を求める。
J(x(k))Δx(k) =
-f(x(k))
5 (ﾌﾟﾛｸﾞﾗﾑ590-650)
x(k+1)
=
x(k)
+Δx(k)
とおく。
6
(ﾌﾟﾛｸﾞﾗﾑ660-670) Δx(k)が収束判定値εより小さいならx(k+1)を方程式f(x)=0の解として計算を終了する。そうでないならx(k+1)を新しい近似値として2に戻る。
BASICによる連立非線形方程式解法プログラム
100 '==連立非線形方程式解法プログラム=====（プログラムの大きさ３．２Ｋ）
110 ' Pseudo Newton's Method
120 ' by ｽﾉｳﾁ ﾊﾙｵ " ﾊﾟｿｺﾝﾆﾖﾙｽｳﾁｹｲｻﾝ ",p.144,ｻｲｴﾝｽｼｬ
210 READ NN
220 DIM ZA(NN, NN), ZB(NN), ZX(NN), ZY(NN), NG(NN)
230 DIM F(NN), ZF1(NN), ZX1(NN), X(NN)
240 '-----------------
280 FOR NI = 1 TO NN: READ ZX1(NI): NEXT
300 NL1 = 100
310 ZE =
.00001 （これは1変数非線型方程式の解法のNewton-Raphson
320 ZD = .001 法をn変数に拡張したものである
330 '---------------- 3080行でJacobi行列を前進差分近似で求めている
340 FOR NL = 1 TO NL1 したがって，初期値がある程度近くないと収束しな
350 FOR NI = 1 TO NN いことに注意。またNewton法同様expに弱い。）
360 X(NI) = ZX1(NI) 近似解を代入
370 NEXT NI
380 '----------------
390 GOSUB 4000
400 FOR NI = 1 TO NN
410 ZF1(NI) = F(NI)
420 NEXT NI
430 '----------------
440 GOSUB 3000
450 '----------------連立一次方程式をGaussの消去法で解く
470 FOR NI = 1 TO NN
480 NG(NI) = NI
490 NEXT NI
500 GOSUB 1000
510 '----------------
520 FOR NI = 1 TO NN
530 NJ = NG(NI): ZB(NI) = -ZF1(NJ)
540 NEXT NI
550 GOSUB 2000
560 '----------------
580 PRINT "ﾊﾝﾌﾟｸ ｶｲｽｳ "; NL
590 NC = 0
600 FOR NI = 1 TO NN
610 ZW = X(NI)
620 ZX1(NI) = ZX1(NI) + ZW
630
IF ABS(ZW) < ZE OR ABS(ZW) < ZE * ABS(ZX1(NI)) THEN NC = NC + 1
640 PRINT "X("; NI; ")="; ZX1(NI)
650 NEXT NI
660 '----------------
670 IF NC = NN THEN 730
680 NEXT NL
690 '----------------
700 BEEP: PRINT "ﾊﾝﾌﾟｸ="; NL1; "ﾃﾞ ｼｭｳｿｸｼﾅｲ"
710 GOTO 770
720 '----------------
730 PRINT "*****ﾎｳﾃｲｼｷ ﾉ ｶｲ*****"
740 FOR NI = 1 TO NN
750 PRINT " X("; NI; ")= "; ZX1(NI)
760 NEXT NI
770 END
1000 '---------------
1010 '---------------
1020 FOR NK = 1 TO NN - 1
1030 '-------------
1040 ZW = ABS(ZA(NK, NK)): NK1 = NK
1050 FOR NI = NK + 1 TO NN
1060 ZW1 = ABS(ZA(NI, NK))
1070 IF ZW1 < ZW THEN 1090
1080 ZW = ZW1: NK1 = NI
1090 NEXT NI
1100 '-------------
1110 IF NK = NK1 THEN 1160
1120 FOR NJ = 1 TO NN
1130 ZW = ZA(NK, NJ): ZA(NK, NJ) = ZA(NK1, NJ): ZA(NK1,
NJ) = ZW
1140 NEXT NJ
1150 NI1 = NG(NK): NG(NK) = NG(NK1): NG(NK1) = NI1
1160 '-------------
1170 ZW = ZA(NK, NK)
1180 FOR NI = NK + 1 TO NN
1190 ZA(NI, NK) = ZA(NI, NK) / ZW
1200 FOR NJ = NK + 1 TO NN
1210 ZA(NI, NJ) = ZA(NI, NJ) - ZA(NI, NK) * ZA(NK, NJ)
1220 NEXT NJ
1230 NEXT NI
1240 NEXT NK
1250 RETURN
2000 '-------SOLVE
2010 '---------------
2020 ZY(1) = ZB(1)
2030 FOR NI = 2 TO NN
2040 ZW = ZB(NI)
2050 FOR NK = 1 TO NI - 1
2060 ZW = ZW - ZA(NI, NK) * ZY(NK)
2070 NEXT NK
2080 ZY(NI) = ZW
2090 NEXT NI
2100 '----------------
2110 X(NN) = ZY(NN) / ZA(NN, NN)
2120 FOR NI = NN - 1 TO 1 STEP -1
2130 ZW = ZY(NI)
2140 FOR NK = NI + 1 TO NN
2150 ZW = ZW - ZA(NI, NK) * X(NK)
2160 NEXT NK
2170 X(NI) = ZW / ZA(NI, NI)
2180 NEXT NI
2190 RETURN
3000 '----------------Jacobi行列の計算
3010 '----------------
3020 FOR NJ = 1 TO NN
3030 IF NJ <> 1 THEN X(NJ - 1) = ZX1(NJ - 1)
3040 X(NJ) = X(NJ) + ZD
3050 '--------------
3060 GOSUB 4000
3070 FOR NI = 1 TO NN
3080 ZA(NI, NJ) = (F(NI) - ZF1(NI)) / ZD
3090 NEXT NI
3100 NEXT NJ
3110 RETURN
4000 '**** ﾎｳﾃｲｼｷ 関数値ベクトルf(x)の計算
4005 '
4010 '(Z,N ﾉﾂｸﾍﾝｽｳﾊﾂｶﾜﾅｲ )（F(1)=,F(2)=...
の右辺に
4020 '--ﾐﾁｽｳ ﾉ ｶｽﾞ X() f1(x1,x2,...)=0
4030 DATA 4
f2(x1,x2,...)=0
4040 '--ｼｮｷﾁ X(1),X(2)... .....
4051 DATA 350,12,5,30 fn(x1,x2,...)=0
4060 '--ﾚﾝﾘﾂ ﾎｳﾃｲｼｷ の左辺を書く）
4100 F(1) = X(3) + 3 * X(4) - 90
4110 F(2) =X(3)+X(1)-(90 +X(1)) * (1-0.??)
一回通過転化率↑
4120 F(3) = X(1) - X(3) * X(2) / 1
4130 F(4) = 90*(1-0.??)-X(3)
↑総括収率
4220 RETURN
表計算の機能「ソルバー」による連立非線形方程式の解法
表計算の機能「ソルバー」は変数が複数の連立非線形方程式を解くための道具である。しかし収束させる値はひとつなので使い方に工夫が必要である。普通は「連立方程式の残差二乗和を最小にする」という条件により計算させる。 (pr07.xls)
【例題】（化学工学基礎のリサイクル－パージ問題より）
石炭の水蒸気存在下の部分酸化によって生成した COとH2を洗浄したあと燐酸触媒反応器に送り，メタノール合成反応： CO+2H2→CH3OH
をおこなう。原料ガスはCO,H2が量論比の組成(1:2)で供給される。反応器におけるCOの一回通過転化率は18%である。および総括収率
K=0.95 の条件でプロセス中の各流量を計算する。
【解】
上図のように各流量を未知数ｘ1，ｘ2，ｘ3，ｘ4とする。以下のように物質収支式をつくる。
①系全体の収支 
②反応器まわりり収支 （原料）x( 1-（一回通過転化率）) ＝（未反応物）:
③制約条件１ ： リサイクル流れとパージ流れは組成が同じ:
④制約条件２ 総括収率Ｋ=0.95：
この連立方程式は③の条件で非線形となっている。以下のように表計算により解く。(pr07.xls)
化学工学における図式解法と非線形連立方程式
精留塔（物質収支と平衡による非線形連立方程式）
実際の蒸留塔の濃縮部を模した４段の段塔により蒸留をおこなう．原料は流量V,濃度 zFの蒸気として塔底部に供給される．これを還流比 Rで精留をおこない，流量D,濃度 xDの留出液を得る．この装置は回分式であるため，実際には留出液の抜き出しで zFは変化するが，ここでは zF一定の定常状態にあるとして解析する．
蒸気流量V[mol/s],液流量L[mol/s]は一定と仮定する．全体の物質収支よりV=D+L, および zFV=DxD+Lx4が成立する．還流比(Reflux ratio)RをR=L/Dで定義する．8個の未知数 xD,x1,y2,x2,y3,x3,y4,x4についての方程式は以下の８個である．
①全体の物質収支：
②塔頂から各段の下までの物質収支：
（式(2)～(3)はyn+1とxnの関係を与え，x-y座標上では直線となる．これが操作線である.）
③各段を去る気液が平衡にあるという条件すなわち理論段の仮定より各段で次式が成立する．
（α=3とする）
【実習】以上の連立方程式をゴールシークで解く。図式解法と比較する。(ex_dis4.xls)
多重効用蒸発（物質収支と熱収支による非線形連立方程式）
３重効用蒸発缶（順流）を用いてＴF= 303K,ｘF=0.05重量分率，のある水溶液 Ｆ= 3.0 kg/s を xW = 0.5 重量分率まで濃縮する。各缶の伝熱面積 は等しくＡ [m2] 、および総括伝熱係数は各々，Ｕ1=2600, Ｕ2=2100, Ｕ3=1400 W/(m2･K) であり，蒸発潜熱，熱容量も温度依存性を無視して各々ｒ=2.3×106J/kg,ｃp=4.2×103J /(kg･K) とする。ＴB3=325K,Ｔs= 403Kとして伝熱面積Ａ[m2]，温度ＴB1,ＴB2,各缶の蒸発量Ｖ1[kg/s],Ｖ2,Ｖ3,供給蒸気量ＶS を求めよ。
解）＜全物質収支＞ Ｆ-Ｗ = Ｖ1 +
Ｖ2 +
Ｖ3
＜溶質収支＞ ｘFＦ =
ｘWＷ
上式より， Ｖ1+Ｖ2+Ｖ3 = ( 1 -ｘF/ｘW )Ｆ
(1)
＜蒸発潜熱に関する熱収支＞
第１缶： ｒＶS =
ｒＶ1 +
Ｆｃp(ＴB1-ＴF)
(2)
第２缶： ｒＶ1
= ｒＶ2 + (Ｆ-Ｖ1)ｃp(ＴB2-ＴB1)
(3)
第３缶： ｒＶ2 =
ｒＶ3 + (Ｆ-Ｖ1-Ｖ2)ｃp(ＴB3-ＴB2) (4)
＜各缶の伝熱速度＞
第１缶： Ｕ1Ａ(Ｔs
-ＴB1) =
ｒＶS
(5)
第２缶： Ｕ2Ａ(ＴB1-ＴB2) = ｒＶ1 (6)
第３缶： Ｕ3Ａ(ＴB2-ＴB3) = ｒＶ2
(7)
以上の７つ未知数による連立方程式を解く。すなわち以下の式の残差を零にする。
f(1)=(1-xF/xW)*F-(V1+V2+V3)
f(2)=r*V1+F*cp*(TB1-TF)
f(3)=r*V2+(F-V1)*cp*(TB2-TB1)-r*V1
f(4)=r*V3+(F-V1-V2)*cp*(TB3-TB2)
f(5)=r*Vs-U1*A*(Ts-TB1)
f(6)=r*V1-U2*A*(TB1-TB2)
f(7)=r*V2-U3*A*(TB2-TB3)
