Python 連立方程式の解(y=f(x)の係数)を求める。そして代入して計算する方法

# 19/10/08刷新
座標2点もしくは3点から、関数y=f(x)へフィッティングするための係数をnumpyで連立方程式を計算する方法です。その後、適当な複数のXをその関数へぶちこんで、X,Yデータをcsvファイル出力するまでを行う。

f:id:HK29:20191008172211p:plain

f:id:HK29:20191008172443p:plain

▼本プログラム

行列含む数値計算はnumpyを使用し、数値データだけでなく行列データを操作するにはpandasを使用することで圧倒的に楽になる。

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd

# 関数A：座標2点から一次関数を求める(Y=ax+b)。numpyを使用しない場合
def funcA(Xa, Ya, Xb, Yb):
    a = (Yb - Ya)/(Xb - Xa)
    b = Yb - a * Xb
    print("############ function A")
    print('a={0}, b={1}'.format(a, b))

    return a, b

def exec_funcA(coefficient, inputX):
    Y = coefficient[0] * inputX + coefficient[1]
    print("exec_funcA")
    print(Y)
    
    return Y 

# 関数B：座標2点から一次関数を求める(Y=ax+b)。numpyを使用した場合
def funcB(Xa, Ya, Xb, Yb):
    myXset = np.array([Xa, Xb]) # X座標群をnumpyを利用しての行列格納
    myYset = np.array([Ya, Yb]) # Y座標↑
    coefficient = np.polyfit(myXset, myYset, 1) # 係数a, bを計算してる
    print("############ function B")
    print('a, b =' + str(coefficient))
    Y_func = np.poly1d(coefficient) # y=f(x) 一次元多項式を作成できる
    print(Y_func)

    return Y_func

# 関数C：座標3点から二次関数を求める(Y=ax2+bx+c)。numpyを利用
def funcC(Xa, Ya, Xb, Yb, Xc, Yc):
    myXset = np.array([Xa, Xb, Xc])
    myYset = np.array([Ya, Yb, Yc])
    coefficient = np.polyfit(myXset, myYset, 2) # 係数a, b, cを計算してる
    print("############ function C")
    print('a, b, c =' + str(coefficient))
    Y_func = np.poly1d(coefficient)
    print(Y_func)
    
    return Y_func

def exec_funcBorC(Y_func, inputX):
    Y = Y_func(inputX)
    print("exec_funcBorC")
    print(Y)

    return Y

def writeFunc(outname, myVal):
    tmp_list = ['Y_' + outname + '=', str(myVal), '\n']
    with open("output.txt", "a") as f:
        f.writelines(tmp_list)
        

# 得られた関数Y=f(x) にxを代入して、結果Yを取得する。
def calc(filename, myfunc, X_list):
    print('substitute X for Y in the equation')
    print("set X")
    # numpy形式(array)に変換
    X_np = np.array(X_list)
    print(X_np)
    # Yを計算
    Y_np = myfunc(X_np)  # 複数のXに対するそれぞれのYを一気に計算できる。
    print(Y_np)
    
    #XとYを結合。1次元の縦結合せずに2次元行列にするために、一旦2次元に変換する
    X_np2 = X_np[None,:].T # 転置するため.Tを付ける
    Y_np2 = Y_np[None,:].T
    # numpy形式の2次元配列の結合
    XY_np = np.concatenate([X_np2, Y_np2], axis=1)
    print("XY_numpy")
    print(XY_np)
    # リスト(普通の配列形式)へ変換
    XY_list = XY_np.tolist() # tolistメソッドを利用
    print("XY_list")
    print(XY_list)
    # pandas形式(DataFrame)へ変換
    df = pd.DataFrame(data=XY_list, columns=['x', 'y'])
    print("XY_pandas")
    print(df)
    # csvで保存する
    df.to_csv(filename + "_output2.csv", header=True, index=False, encoding="utf-8")
    

if __name__ == '__main__': #このファイルを本体として実行した場合、mainが実行される。
    Xa = 15.8
    Ya = 0
    Xb = 23.4
    Yb = -1.2e18
    Xc = 25.8
    Yc = -1.7e18
    inputX = 27
    X_list = [-2, 0, 1.58, 2.1, 2.34, 3, 7, 14.2, 18, 22, 27, 30]
    
    # 関数Aの作成と結果抽出
    coefficient = funcA(Xa, Ya, Xb, Yb)
    outputY = exec_funcA(coefficient, inputX)
    writeFunc("funcA", outputY)
    
    # 関数Bの作成と結果抽出
    myfuncB = funcB(Xa, Ya, Xb, Yb)
    outputY = exec_funcBorC(myfuncB, inputX)
    writeFunc("funcB", outputY)
    
    # 関数Cの作成と結果抽出
    myfuncC = funcC(Xa, Ya, Xb, Yb, Xc, Yc)
    outputY = exec_funcBorC(myfuncC, inputX)
    writeFunc("funcC", outputY)
    
    # 作成した関数にXを代入してYを計算
    myXY = calc("myfuncB", myfuncB, X_list)
    myXY = calc("myfuncC", myfuncC, X_list)