Python で ode45()関数を模倣する

Aditya Raj 2023年10月10日
Python で ode45()関数を模倣する

常微分方程式は、多くの科学的問題を解決するために MatLab で使用されます。ode45() は、MatLab で微分方程式を解くために使用されます。

この記事では、Python で ode45() 関数を模倣する方法を説明します。

Python で ode45() 関数を模倣する

Python で ode45() 関数を模倣するために、scipy モジュールで定義された solve_ivp() メソッドを使用できます。solve_ivp() メソッドは、常微分方程式(ODE)のシステムを統合します。

  • solve_ivp() メソッドは、最初の入力引数として関数を取ります。input 引数で指定された関数は、微分方程式の係数を含む配列を返す必要があります。
  • -2 番目の入力引数では、solve_ivp() メソッドは 2つの数値を含むタプルまたはリストを取ります。値は積分の間隔を表し、タプルの最初の値は間隔の開始を表し、タプルの 2 番目の値は間隔の最大値を表します。
  • 3 番目の入力引数では、solve_ivp() メソッドは初期値を表す配列を取ります。
  • 実行後、solve_ivp() メソッドはさまざまな属性を持つ bunch オブジェクトを返します。
    1.t 属性には、時点を含む numpy 配列が含まれています。
    2.y 属性には、t の値と時点を持つ numpy 配列が含まれています。
    3.sol 属性には、微分方程式の解を含む Odesolution オブジェクトが含まれています。solve_ivp() メソッドで dense_output 引数が false に設定されている場合、sol 属性には None が含まれます。

これをよりよく理解するには、次の例を参照してください。

from scipy.integrate import solve_ivp


def function(t, y):
    return 2 * y


interval = [0, 10]
initial_values = [10, 15, 25]
solution = solve_ivp(function, interval, initial_values)
print("Time:", solution.t)
print("Y:", solution.y)

出力:

Time: [ 0.          0.07578687  0.56581063  1.18741382  1.85887096  2.55035821
  3.25007544  3.95320486  4.65775424  5.36289544  6.06828346  6.77377445
  7.47930839  8.18486026  8.89041961  9.59598208 10.        ]
Y: [[1.00000000e+01 1.16366412e+01 3.10073783e+01 1.07492109e+02
  4.11689241e+02 1.64114780e+03 6.65071446e+03 2.71362627e+04
  1.11036049e+05 4.54874443e+05 1.86437495e+06 7.64300835e+06
  3.13352156e+07 1.28474398e+08 5.26752964e+08 2.15973314e+09
  4.84541488e+09]
 [1.50000000e+01 1.74549617e+01 4.65110674e+01 1.61238163e+02
  6.17533861e+02 2.46172171e+03 9.97607169e+03 4.07043941e+04
  1.66554074e+05 6.82311665e+05 2.79656243e+06 1.14645125e+07
  4.70028233e+07 1.92711598e+08 7.90129446e+08 3.23959970e+09
  7.26812231e+09]
 [2.50000000e+01 2.90916029e+01 7.75184457e+01 2.68730272e+02
  1.02922310e+03 4.10286951e+03 1.66267862e+04 6.78406569e+04
  2.77590123e+05 1.13718611e+06 4.66093739e+06 1.91075209e+07
  7.83380389e+07 3.21185996e+08 1.31688241e+09 5.39933284e+09
  1.21135372e+10]]

上記の例では、最初に、入力引数として ty を取り、y に基づいて値を返す function という名前の関数を定義しました。

次に、変数 intervalinitial_values をそれぞれ使用して、ODE の間隔と初期値を定義しました。functioninterval、および initial_valuessolve_ivp() 関数への入力引数として渡し、最後に、変数ソリューションで出力を取得します。

出力では、時間値が 0 から 10 の間隔全体に広がっていることがわかります。同様に、出力には、各時間値に対応する y 値が含まれています。

ソリューションの属性 t で時点を明示的に指定することもできます。このために、以下に示すように、y 値が必要な目的の時間値を含む配列を solve_ivp() メソッドの t_eval 引数に渡す必要があります。

from scipy.integrate import solve_ivp


def function(t, y):
    return 2 * y


interval = [0, 10]
initial_values = [10, 15, 25]
time_values = [1, 2, 3, 6, 7, 8]
solution = solve_ivp(function, interval, initial_values, t_eval=time_values)
print("Time:", solution.t)
print("Y:", solution.y)

出力:

Time: [1 2 3 6 7 8]
Y: [[7.38683416e+01 5.46053271e+02 4.03089733e+03 1.62618365e+06
  1.20160156e+07 8.87210156e+07]
 [1.10802512e+02 8.19079906e+02 6.04634600e+03 2.43927547e+06
  1.80240234e+07 1.33081523e+08]
 [1.84670854e+02 1.36513318e+03 1.00772433e+04 4.06545912e+06
  3.00400390e+07 2.21802539e+08]]

時間値には、t_eval パラメーターへの入力引数として渡される値のみが含まれていることがわかります。同様に、属性 y には、指定された t 値のみの値が含まれます。

このアプローチは、間隔内の特定のポイントの値を取得するのに役立ちます。

著者: Aditya Raj
Aditya Raj avatar Aditya Raj avatar

Aditya Raj is a highly skilled technical professional with a background in IT and business, holding an Integrated B.Tech (IT) and MBA (IT) from the Indian Institute of Information Technology Allahabad. With a solid foundation in data analytics, programming languages (C, Java, Python), and software environments, Aditya has excelled in various roles. He has significant experience as a Technical Content Writer for Python on multiple platforms and has interned in data analytics at Apollo Clinics. His projects demonstrate a keen interest in cutting-edge technology and problem-solving, showcasing his proficiency in areas like data mining and software development. Aditya's achievements include securing a top position in a project demonstration competition and gaining certifications in Python, SQL, and digital marketing fundamentals.

GitHub

関連記事 - Python Function