csvファイルを用いて関数を定義したいです

Question

#####要旨
Pythonでローレンツ方程式のような常微分方程式をscipy.integrate.odeintなどを用いて解く場合に、
関数の定義をcsvファイルなどを利用して簡便に行いたいです。

#####詳細1
・scipy.integrate.odeintの使い方まとめ：
https://qiita.com/Sunset_Yuhi/items/d938718ad277eeab746c
を参考にすると、
ローレンツ方程式
【式１】

を

引用テキスト

def func_lorenz(var, t, p, r, b):
    dxdt = -p*var[0] +p*var[1]
    dydt = -var[0]*var[2] +r*var[0] -var[1]
    dzdt = var[0]*var[1] -b*var[2]

    return [dxdt, dydt, dzdt]

といった具合で直接記述しています。
この方法だと、方程式の項が多い場合は入力に時間がかかってしまい、不便を感じています。
そこで関数の定義をcsvファイルなどを利用して簡便に行いたいと考えたのですが、方法がわからず、この度質問いたしました。

例えば、ローレンツ方程式を
【式２】

のように行列表現にし、
Aの部分とX'の部分をcsvファイルか何かで用意して、
pd.read_csv()などで読み込ませれば、何かできますでしょうか。

引用しました
＞ dydt = -var[0]var[2] +rvar[0] -var[1]
の
＞var[0]*var[2]
のような部分をcsvファイルなどを用いて、どのように記述すればよいか見当がつかず、
教授いただけると幸いです。

まずはパラメーターp,r,bを与えた、
例えば
【表1】

のようなcsvファイルで行列Aを定義することを考えています。
対応する変数X'は、どのようにcsvファイルで用意すればよいかわからず、、
例えば
【表2】

のように2列に分けて記述し、A列×B列という処理をpython側で行えば何かできますでしょうか。

舌足らずかもしれず大変恐縮ですが、よろしくお願いいたします．

#####詳細2
詳細1に記述した、
行列Aの一部のパラメーターを時間変動させる場合の関数の定義方法を知りたいです。
例えば
パラメーターpを
【表3】

といったようにcsvファイルで用意しておき、
行列Aを
【表4】

のようにcsvで用意しておき、
python側で図3と図4のpを関連付けるなど、できるのでしょうか。

何か良い方法がありましたら、ご教授いただけると幸いです。
よろしくお願いいたします。

#####追記
toast-uz様の回答を参考に、【詳細２】に対応するコードを考えてみました。
ベストアンサーはtoast-uz様につけさせていただいたので、こちらに追記することにしました。

詳細１の式１のAをmatとして
【表５】

のように行列で定義しました。

さらに、パラメーターpの時間変動分を
【表６】

のようにしました。
表６の列0はt、列1はpです。
行数　＝　計算終了時間/計算間隔にしています。

関数部分は

def func_lorenz(var, t, mat, p_list):
    x, y, z = tuple(var)
    mat[0,0] = - p_list[p_list[:,0] == round(t,2), 1]
    mat[0,1] = p_list[p_list[:,0] == round(t,2), 1]
    dXdt = [0,0,0]
    for i in range(len(dXdt)):
        dXdt[i] = mat[i,0]*x + mat[i,1]*y + mat[i,2]*z + mat[i,3]*x*y +  mat[i,4]*x*z
    return dXdt

のように定義しました。
パラメーターpに対応するmat[0,0]とmat[0,1]を」関数内で書き換えるようにしました。
また、連立数が多い方程式を解く場合も想定し、for分を使って連立方程式を記述しました。
必要な基底関数のみ直接記述し、対応する係数をcsvなどで用意する方法をとりました。

一応計算はできましたが、効率の悪いコードになっていると思っております。
もし、改善するところがありましたらコメントいただけると大変助かります。

Comments

[toast-uz]

関数の引数を7次元ベクトルα=(x, y, x, t, p, r, b)T、戻り値を3次元ベクトルβ=(dxdt, dydt, dzdt)Tとしたとき（Tは転置）に、戻り値を引数から求める計算式を作れますか？例示されたように行列内に未知数を入れたり、xyなどと掛け算の項を入れてはいけません。未知数は必ず上の7次元ベクトルだけを使って（何回も使って良いです）、数値だけの行列とあわせた計算により、戻り値を求められるようにしてください。それができるなら、実現可能だと思います。

なんとなく、β = B×α'×A×α'T　（Tは転置、Aは数値のみの8×8行列、Bは数値のみの3×8行列、α'はαの末尾に1を追記した8次元ベクトル）という形に書けそうではあります。

これができるなら、汎用的な、引数n次元、戻り値m次元、A=(n+1)×(n+1)行列とB=m*(N+1)行列をcsvで与えると、引数の任意の2次以下の関数として、戻り値を表すことができます。

ほぼ回答になってしまいましたので、こちらで回答を書きます。

[toast-uz]

ちょっと行列の考え方が違っていたので回答蘭で修正します。

[tmyuki_0716]

ありがとうございます。
回答を拝見いたしました。

[toast-uz]

追記部分、dXdt = [0,0,0] 以下returnまでを

dXdt = mat @ np.array([[x, y, z, x*y, x*z]]).T
return dXdt.ravel().tolist()

と表現することができます。numpy利用時は、forループをもっとうまく書けないか工夫しましょう。
@は行列の積、Tは転置です。こんな風に行列計算できると、カッコ良いですね。

[tmyuki_0716]

ありがとうございます。
こんなに簡単にかけるのですね！
勉強になりました。

Answer 1

質問者様と議論したところ、最終的な拡張を考えると、csv等に定義を外だしするのは複雑化してしまいそうです。むしろ、既存の直接記述を簡単にする方が早そうだ、となりました。よって、その方針の回答を記載します。

以下のように引数varの配列をタプルに変換して個々の変数として表現してあげると、直感的でわかりやすい記述になります。

def func_lorenz(var, t, p, r, b):
    x, y, z = tuple(var)
    dxdt = -p*x + p*y
    dydt = -x*z + r*x - y
    dzdt = x*y - b*z
    return [dxdt, dydt, dzdt]

Comments

[tmyuki_0716]

ありがとうございます。
大変助かりました。
また、勉強になりました。
最後に、
toast-uz様の回答を参考に、詳細２に対応するコードを考え、
質問本文に追記しました。
何か改善案がありましたら、コメントいただけると幸甚です。
様々ありがとうございました。

Answer 2

コード化はせずに回答方針のみ書きますので、コード化はやってみてください。

まず、質問者様が例示した表は、いずれも未知数を含むものですので、「未知数を解釈して計算する」という大きなハードルが発生します。sympyを使う方法もありそうですが、sympyでの数式定義は複雑であり、質問者様の要件（簡単に数式を定義したい）に対しては本末転倒になってしいます。よって**、「未知数をうまく分離して任意の数式を表現する」ことが必要です**。これができれば、質問者様の目的は簡単な実装に帰着します。

以下、「未知数をうまく分離して任意の数式を表現する」ための方式を考えてみましょう。

関数の引数をm次元の列ベクトルα、戻り値をn次元の列ベクトルβとします。
質問の場合は、引数は7次元ベクトルα=(x, y, z, t, p, r, b)T、戻り値は3次元ベクトルβ=(dxdt, dydt, dzdt)Tです。この場合、dxdtは1つの変数とみなしxの微分がーとかは考えません。Tは転置ですので、それぞれ列ベクトルを表現しています。

α'を、αの最後に1を付加したm+1次元の列ベクトルとします。

そのとき、数値からなる行列Ai(m+1行m+1列) （iは添字で、0<=i<n）をうまくとれば、任意の、「引数の2次以下の式」は

α'T×Ai×α' （Tは転置、iは添字で0<=i<n）

で表現できます。αではなくα'を使うことがポイントで、これにより、2次未満の数式項を表現することができます。ちなみに、この時、Aは上三角行列として表現でき、各要素の位置に対応した未知数（もしくは1）の掛け算の係数を表します。ですので、数式をもとに行列Aは自明に求めることが可能であり、実用上の使い勝手の良い形であると言えます。

よって、β=(β1, β2, ・・・)とすると、

βi = α'T×Ai×α' （Tは転置、iは添字で0<=i<n）

となり、「未知数をうまく分離して任意の数式を表現する」ことが達成できました。

Pythonは辞書型やリスト型を活用することで、引数も戻り値も任意の次元を表すことができます。よって、行列Aiをcsv等で入力することで、簡単に上記計算ができます。質問者様の目的は達成です。

Comments

[tmyuki_0716]

ご教授いただきありがとうございます。
丁寧に説明いただいたおかげで、理解することができました。
コード化については不慣れなので、色々調べながらトライしてみます。
できましたら、改めてこの場でコメントいたします。

詳細２の件は、行列Aの中で値が変動する要素を引数tの値に従って関数内で書き換える方法を思いつきました。これについてもコード化にトライしてみます。

また、最終的に
dx/dt = a x^4 + b(x-y)^0.8+・・・
のように4乗程度の高い指数を扱ったり、
変数同士を引いた値に1未満の指数を与えたり、
を実現することを目指しておりますが、
回答者様の方法を応用することで実現可能でしょうか。
実現可否だけでも教授いただけると幸甚です。

[toast-uz]

不可能です。

まず4乗であっても多項式で表せる範囲であれば、理論的には行列Aを2次元から多次元に拡張していけば表現可能です。しかし、その行列はとてもスパースになり、無駄が多く、直感的にもわかりにくいです。そのため、「直接記述より簡単にしたい」という当初目的の達成ができなくなります。ましてや、新しく例示いただいたb(x-y)^0.8は多項式で表現できません。今後、三角関数や対数関数の項も出でくるのではないでしょうか。

また、sympy含め、式の構文を表現して解釈するような仕組みを入れたとしても、その式を記述するくらいでしたら、最初の直接記述の方が手っ取り早いです。

以上に対して汎用的なアプローチはなく、直接記述することが望ましいです。

[tmyuki_0716]

ありがとうございます。勉強になりました。
確かにスパースなので、
直接記述で検討しようと思います。

[toast-uz]

根本的な疑問なのですが、var[i]という表現にしているから「面倒感」があるものと思います。
以下の表現だとかなり簡素ですが、これではダメなのでしょうか？
def func_lorenz(x, y, z, t, p, r, b):
dxdt = -p*x + p*y
dydt = -x*z + r*x - y
dzdt = x*y - b*z
return [dxdt, dydt, dzdt]

もしくはvarのままでも、以下のような記述が可能です。
def func_lorenz(var, t, p, r, b):
x, y, z = tuple(var)
dxdt = -p*x + p*y
dydt = -x*z + r*x - y
dzdt = x*y - b*z
return [dxdt, dydt, dzdt]

[tmyuki_0716]

早速コメント頂きありがとうございます。
＞x, y, z = tuple(var)
という方法もあるのですね。

勉強になります。
確かに変数の名前が長すぎるのが原因かもしれません。
教示いただいた表現方法でトライしてみます。

[toast-uz]

それでは、いちおう、こちらの方法も回答として別に記載しておきます。