数学建模 |
您所在的位置：网站首页 › scipy解方程 › 数学建模 |
|
introduction：
python对于常微分方程的数值求解是基于一阶方程进行的，高阶微分方程必须化成一阶方程组，通常采用龙格-库塔方法. scipy.integrate模块的odeint模块的odeint函数求常微分方程的数值解，其基本调用格式为： sol=odeint(func,y0,t) func是定义微分方程的函数或匿名函数 y0是初始条件的序列 t是一个自变量取值的序列(t的第一个元素一定必须是初始时刻) 返回值sol是对应于序列t中元素的数值解，如果微分方程组中有n个函数，返回值会是n列的矩阵，第i(i=1,2···,n)列对应于第i个函数的数值解. 其使用的方法和栗子我们可以查看官方给的栗子： Examples -------- The second order differential equation for the angle `theta` of a pendulum acted on by gravity with friction can be written:: theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0 where `b` and `c` are positive constants, and a prime (') denotes a derivative. To solve this equation with `odeint`, we must first convert it to a system of first order equations. By defining the angular velocity ``omega(t) = theta'(t)``, we obtain the system:: theta'(t) = omega(t) omega'(t) = -b*omega(t) - c*sin(theta(t)) Let `y` be the vector [`theta`, `omega`]. We implement this system in python as: >>> def pend(y, t, b, c): ... theta, omega = y ... dydt = [omega, -b*omega - c*np.sin(theta)] ... return dydt ... We assume the constants are `b` = 0.25 and `c` = 5.0: >>> b = 0.25 >>> c = 5.0 For initial conditions, we assume the pendulum is nearly vertical with `theta(0)` = `pi` - 0.1, and is initially at rest, so `omega(0)` = 0. Then the vector of initial conditions is >>> y0 = [np.pi - 0.1, 0.0] We will generate a solution at 101 evenly spaced samples in the interval 0 > t = np.linspace(0, 10, 101) Call `odeint` to generate the solution. To pass the parameters `b` and `c` to `pend`, we give them to `odeint` using the `args` argument. >>> from scipy.integrate import odeint >>> sol = odeint(pend, y0, t, args=(b, c)) The solution is an array with shape (101, 2). The first column is `theta(t)`, and the second is `omega(t)`. The following code plots both components. >>> 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() """ 例题1：给定方程，求在1 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |