と、格好つけてみましたけれども、大袈裟な物言いの割には、小ネタであったりします。
以前、fx-5800P に、内蔵公式集が 128本くらい入っていると知り、更に、ユーザーによる公式も登録可能、との事で、こうした公式集をダラリと使っていたら、そうした公式に少しは理解が及ぶ所で、初学者の助けになるんじゃないのか、などとボンヤリ思う所でした。
しかし、公式集をみると、例えば、multiple equation solver の様なメカがないため、公式内の変数に対して、式を変形したものが重複して登録されております。
理想を言うならば、公式を一本、入れておいて、Solver が指定した変数に対して数値解を求めてくれる、のが欲しいなぁ、というヘタレぶり。
fx-5800P の様な、公式内蔵のプロ電は、目下、successor が存在せず、というのが実に押井、と思う次第。
fx-CG50 には、Solver Appsがあり、方程式を登録する事で、所望の変数について数値解を計算する事は可能であります。
そこで、Solver Appsに、内蔵公式 selector の様なメカがあると便利かなァ、と思うのありますが、fx-CG50にはもちろんの事、fx-CG100 にも、そうした便利メカは用意されておりません。そりゃそうだ、学校教育向けでしたから、内蔵公式集の様な「あんちょこ」(「安直」の訛りが語源だそうです)は無いのが当たり前でありました。
斯様に、プロ電が「プロツール」の座位から「学校教育向け」に移って行く今日、プロツールの動きとしては、fx-FD10 が最後の輝きだったのか、などと思う所であります。
fx-FD10 が death-continue に至った今日、fx-5800Pが、プロツールの孤塁を守っている、という認識で良いのかも知れません。
プロツールのsuccessor 登場、という局面は期待できない今日、代替製品として、fx-CG50/fx-CG100 の upython に活路を見出すという試みは「密かなる需要、有用性が在りはしまいか ? 」という淡い期待から思いつくのであります。
幸いにして、fx-CG50/CG100 には、スクリプトを沢山抱えるだけのストレージ容量があります。
内蔵公式集として導入されていなくとも、ユーザーが必要な公式を登録、メニュー形式なりで計算式を選択して計算するといった作業を upython scriptとして記録しておけば、十分に活用できそうです。
実際、fx-5800Pに土木計算ライブラリを書き込んだものを、「専用機」として流通している商売があるとの話。
で、ここまでが長すぎる"前フリ"でありました。以降が本題。
fx-CG50/CG100 の Solver Apps, 結構便利なのですが、「方程式セレクタの様なメカが欲しい」という点もさることながら、f(x)=0 の形式の方程式しか解けない、そこが"押井"かなぁ ... と。
で、本編では、upython の計算能を用いて、2元連立方程式の数値解 Solver を模索しようという試みであります。
2方程式の数値解を求めるには、Newton-Rafson 法が利用できるそうで、ここでは簡易版のスクリプトを提示し、読者の皆様 (そんなに居らんのよなァ)に一考願おうという次第です。
つぎのコード「slv2eq.py」を用意します。
```
# title : slv2eq.py - Newton-Rafson Numeric solver
# begin : 2025-08-28 23:25:08
# mat[2][3] for equation difference matrix
mat = [ [0 for i in range(3)] for i in range(2) ]
eps = 1e-9 # machine epsilon
# pseudo-diff df/dx
def dfx(f, x, y) :
#print('debug ... eps ', eps)
return (f(x+eps, y)-f(x, y))/eps
# pseudo-diff df/dy
def dfy(f, x, y) :
#print('debug ... eps ', eps)
return (f(x, y+eps)-f(x, y))/eps
# newton-rafson solver
def solvre(f, g, x0, y0) :
x, y = x0, y0
for i in range(20) :
print('.', end='')
#report(x, y)
#
mat[0][0] = dfx(f, x, y)
mat[0][1] = dfy(f, x, y)
mat[0][2] = -f(x, y)
#
mat[1][0] = dfx(g, x, y)
mat[1][1] = dfy(g, x, y)
mat[1][2] = -g(x, y)
# matrix solvre
d = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0]
if d == 0 :
d = eps
tx = ( mat[1][1]*mat[0][2]-mat[0][1]*mat[1][2])/d
ty = (-mat[1][0]*mat[0][2]+mat[0][0]*mat[1][2])/d
mat[0][2] = tx
mat[1][2] = ty
x += mat[0][2]
y += mat[1][2]
print(' END')
print('x=',x)
print('y=',y)
print('eq1()=',f(x, y))
print('eq2()=',g(x, y))
return x, y
def report(x, y) :
print('x = {0:.6f} '.format(x), end='')
print('y = {0:.6f} '.format(y))
def inpnum(str):
lin=input(str).split(',')
if len(lin)>1:
num=complex(float(lin[0]),float(lin[1]))
else:
num=float(lin[0])
return num
```
これを import すると、つぎの2つの関数が利用できます。
- slv2eq.inpnum()
数値入力サービス - slv2eq.solvre()
簡易 Newton-Rafson 数値 Solver
で、使い方としては、
- 3*x+2*y-3 = 0
- 3*x-2*y-8 = 0
の様な方程式を解くため、つぎのスクリプトを実行します。
```
from slv2eq import *
# equation #1
def f(x, y) :
return 3*x+2*y-3
# equation #2
def g(x, y) :
return 3*x-2*y-8
# main part
while True :
x=inpnum('x0= ')
y=inpnum('y0= ')
x, y = solvre(f, g, x,y)
```
のようにします。
- def で方程式 = 0 の形式で関数定義をして、
- slv2eq.inpnum() で 数値解探索を始める初期値入力、
- slv2eq.solvre(f, g, x, y) で Newton-Rafson を動かす
のですが、2, 3 の部分がループになっているのは、初期値の入力が簡単に行え、何度もリトライする為であります。
出力例
```
x0= 1
y0= 1
.................... END
x= 1.833333333333334
y= -1.25
eq1()= 0.0
eq2()= 0.0
```
解の探索を始める初期値 x0，y0 を入力すると、20回ほどの計算の後、数値解を吐き出します。
再び、初期値入力になりますが、この例題では、他に答えがないので、[AC」キーを押して終了します。お疲れ様。
上記の例は、線形2元連立方程式のため、fx-CG50 の Solver Appsでも解けるものでしたが、このスクリプトの威力は、つぎの様な式でも数値解を得られる所です。
この方程式の場合、数値解が複素数になってしまうので、素直に解けません。
所が、Newton-Rafson 法では、探索を開始する初期値に複素数を用いると、解けてしまうのですネ。コリャーエエ !!
slv2eq.inpnum() では、実数と複素数の入力に対応しております。
複素数の入力をするには、実部のあと「,」を入れ、続けて虚部を入力します。
例えば、初期値に 1-2j を入力するには、1,-2 という入力をします。
実行例
```
x0= 1
y0= 1
.................... END
x= -1.5216944570145454
y= -0.6115683445992836
eq1()= 11.31555402050879
eq2()= -0.3419836755275885
x0= 1,1
y0= 1
.................... END
x= 3j
y= (-0.8761088075138545-0.5215379877427583j)
eq1()= 0j
eq2()= 0j
```
他にも、色々な式が計算できます。しかし、数値計算で、内部で「数値差分」を使うなどしている所から、ピッタリの数値を得られるとは限りません。飽くまでも、正確な値に"にじり寄る"手段として、活用して戴けると幸いです。
何だか、小ネタと書いた割に、長くなってしまい、申し訳ない。
【お・ま・け】
HP50G Users guide の Page 7-5 (242) に、Manning 方程式の計算例があります。
これを slv2eq.py の応用例として、書いてみました。コードを示しておきます。
```
#
# title : manning.py
# begin : 2024-05-13 14:19:41 , 2022-05-11 16:17
# note : manning formulae equation solving test.
# : ex. cu/n*(y*(b+m*y))**(5/3)/(b+2*y*sqrt(m*m+1))**(2/3)*sqrt(s0) - q = 0
# : y+(q/(y*(b+m*y)))**2/(2*gv) - h0 = 0
# : ans. q = 20.6614376636849357
# : y = 4.9936961327530183
#
import math
import slv2eq
## variables
# constants
h0 = 5.0 # represents the energy head (m, or ft) available for a flow at the entrance to a channel
b = 1.5 # bottom width (m or ft) of channel
m = 1.0 # the side slope (1V:mH) of the cross section.
n = 0.012 # Manning s coefficient, a measure of the channel surface roughness (e.g., for concrete, n = 0.012)
s0 = 0.00001 # slope of the channel bed expressed as a decimal fraction
gv = 32.2 # gravitic accel coefficient
cu = 1.486 # constant for UNIT convinient (ft, yards)
# equation #1
def f(q, y) :
return cu/n * (y*(b+m*y))**(5/3)/(b+2*y*math.sqrt(m*m+1))**(2/3) * math.sqrt(s0) - q
# equation #2
def g(q, y) :
return y + (q/(y*(b+m*y)))**2/(2*gv) - h0
# main part
x, y = slv2eq.solvre(f, g, 2, 1)
```
