MathematicaのNDSolveをPythonで実現したい

前提・実現したいこと

Mathematicaで解いた微分方程式のプログラムを、Pythonコードに書き換えることを目標としています。

f[t]とg[t]を未知関数とした、以下のような連立微分方程式が成立しています。

f''[t] = -(4ft)/4f[t]^2+4g[t]^2+1

g''[t] = (4gt)/4f[t]^2+4g[t]^2

発生している問題

2階微分方程式を1階微分方程式に変換するために、1次導関数を以下のように新しい未知関数(f_1,f_2,g_1,g_2)として定義し、階数下げをしてscipy.integrate.odeintを利用しようと考えています。

f_1=f[t]
f_2=df[t]/dt
g_1=g[t]
g_2=dg[t]/dt

該当のソースコード

Mathematicaの該当のソースコードは以下の通りです。deqは最初に記述したf''[t]とg''[t]を{f''[t],g''[t]}として方程式を立式しています。

Mathematica
1sol = NDSolve[{deq, {f[0], g[0]} == {0, 0}, {f[1], g[1]} == {1, 1}}, {f, g}, {t, 0, 1}]
2ParametricPlot[{f[t], g[t]} /. sol[[1]], {t, 0, 1}]

試したこと

scipyのodeintを利用して数値解を求めるということ、そして1階常微分方程式におけるodeintの立式方法など基本的な部分は理解することができました。しかし、連立かつ変数が4つ(f_1, f_2, g_1, g_2)の微分方程式の場合、どのように関数や初期値（(f[0],g[0])=(0,0)と(f[1],g[1])=(1,1)）を設定してodeintを利用すればいいか分かりません。

最終的には、matplotlibを使ってグラフを描画したいです。

何かアドバイス等ございましたら教えていただきたいです。よろしくお願いいたします。

補足情報

開発環境は、jupyter notebookを利用しています。

etherbeg

2021/06/14 15:11

scipyのodeintを使うために数値解を求めるのではなく、scipyのodeintを使って数値解を求めるのではないですか。Mathematicaのプログラムの置き換えと考えるのではなく、イチからPythonで書くと考えるのが良いと思います。 "scipy.integrate.odeint" でGoogle検索をすると、scipy.integrate.odeintを使ったコードを載せているサイトがいくつもヒットします。それらの多くはmatplotlibを使ったグラフ描画も行なっています。それらのサイトをいくつか開いていって、そこに書かれている解説やコードをじっくり読めば、どのようにコードを書けば良いかのヒントが得られるのではないでしょうか。

komugi_57

2021/06/14 20:31

丁寧にアドバイスをしていただきありがとうございます。scipyのodeintに対して、誤った認識をしておりました。質問内で数値解について述べている部分は修正させていただきます。また、scipy.integrate.odeintのコードを載せているサイトを参考にしながら再考してみようと思います。ご回答いただきありがとうございました。