شبیه سازی سیستم های دینامیکی، بررسی و کنترل آنها بدون حل معادله سیستمی آن بر اساس زمان و پارامترهای وابسته امکان ناپذیر است. محیطهای محاسبات عددی متعددی برای انجام اینکار مثل متلب موجود می باشند که امکانات متنوعی را فراهم آورده اند. زبان برنامه نویسی پایتون با بهره گیری از کتابخانه محاسباتی NumPy و SciPy به ابزار مناسبی برای محاسبات علمی تبدیل شده است. در این نوشته به دنبال حل عددی دسته معادلات دیفرانسیلی با SciPy هستیم.
بسم الله الرحمن الرحیم
برای حل عددی دسته معادلات دیفرانسیلی معمولی با استفاده از کتابخانه SciPy می توانیم از تابع odeint از ماژول integrate از کتابخانه scipy استفاده کنیم. برای استفاده از این تابع ابتدا بایستی معادله دیفرانسیلی مورد نظرمان را به شکل استاندارد در بیاوریم. یعنی تساوی که یک سمت آن دیفرانسیل یعنی dy/dt و سمت دیگر مقادیر آن بر حسب مقادیر متغیرهای اصلی در زمان مشخص t و زمان t باشد. این شکل می تواند در حالت کلی برداری باشد که در این صورت ما مشغول حل دسته معادلات دیفرانسیلی معمولی می باشیم. مثال زیر نحوه تعریف دو معادله دیفرانسیلی را برای مثال نشان می دهد:
def f(x, t):
x1 = x[0]
x2 = x[1]
x1_d = x2 ** 3. + 1.
x2_d = 1. - x1 ** 3.
return [x1_d, x2_d]
که بیانگر معادله دیفرانسیل زیر است:
dx1/dt = x2 ^ 3 + 1 dx2/dt = 1 - x1 ^ 3
شکل تعریف این تابع حالت استانداردی است که برای استفاده از تابع odeint باید رعایت شود. علاوه بر دو متغیر x و t می توان پارامترهای دیگری را نیز بعنوان ورودی تابع در نظر گرفت ولی لزوما دو پارامتره ابتدایی باید حالت سیستم و زمان باشند. حال ورودی های مورد نیاز تابع odeint را ببینیم.
- func که ورودی اول تابع می باشد همان تابعی است که ما تعریف کرده ایم یعنی معادله دیفرانسیلی است که می خواهیم حلش کنیم.
- y0 دومین ورودی می باشد که حالت اولیه معادله دیفرانسیلی است.
- t ورودی سوم و اندیسهای زمانی است که معادله دیفرانسیل باید در آنها به شکل عددی حل شوند. (دقت بیشتر در این پارامتر باعث دقت بیشتر در حل معادله دیفرانسیل می شود که نحوه تعریف آن را در ادامه با یک مثال خواهیم دید.)
این سه پارامتر ورودی های اساسی تابع odeint می باشند که باید دقیق و مطمئن مقداردهی گردند. پارامترهای دیگری نیز وجود دارند که برای تنظیم دقت خروجی مفید هستند که در مجال دیگر به آنها پرداخته خواهد شد. از بین این پارامترها مهمترینشان چندتایی args می باشد. در صورتی که تابع نمایندگی کننده معادله دیفرانسیل شما غیر از حالت و زمان پارامتر ورودی دیگری نیز داشته باشد، باید آن پارامترهای دیگر به شکل یک چندتایی جمع شده و هنگام فراخوانی تابع odeint بعنوان پارامتر args ورودی داده شود.
حال به مثال زیر توجه نمایید:
import numpy as np
from scipy.integrate import odeint
def pend(y, t, b, c):
theta, omega = y
dydt = [omega, -b*omega - c*np.sin(theta)]
return dydt
b = 0.25
c = 5.0
y0 = [np.pi - 0.1, 0.0]
t = np.linspace(0, 10, 101)
sol = odeint(pend, y0, t, args=(b, c))
import matplotlib.pyplot as plt
plt.plot(t, sol[:, 0], 'b', label='theta(t)')
plt.plot(t, sol[:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
این کد با حل معادله پاندول معکوس برای پارامترهای مشخص خروجی حالتها را رسم می نماید. پارامتر t بعنوان ورودی تابع odeint به شکل اعدادی گسسته در بازه ای بین صفر و ده که به شکل خطی در بازه پخش شده اند، ساخته شده است. برای این کار از تابع linspace استفاده شده است که بازه ۰ تا ۱۰ را به ۱۰۱ قسمت یعنی با دقت یک دهم تقسیم می نماید. نمودار خروجی این کدها به شکل زیر می باشد:
