数学建模-用Sympy库求解一些微分方程
例题1:解下列微分方程的特解
from sympy.abc import x #因为sympy的计算都是需要符号变量声明的 例如x=symbols('x'); 使用abc库可以将任何拉丁和希腊字母导出为符号型的
from sympy import diff, dsolve,simplify,Function
y=Function('y') #声明y代表一个函数
eq=diff(y(x),x,2)+2*diff(y(x),x)+2*y #定义方程
con={
y(0):0,diff(y(x),x).subs(x,0): 1} #定义初值条件
y=dsolve(eq,ics=con)
print(simplify(y)) #simplify 化简函数表达式
对于dsolve函数–
def dsolve(eq, func=None, hint="default", simplify=True,
ics= None, xi=None, eta=None, x0=0, n=6, **kwargs):
其参数的含义:eq就是任何常微分形式的等式(甚至可以是微分方程组或者矩阵形式),其表达式都被假设=0
ics就是初始值和边界值的预定义,一般形式如下: ``{f(x0): x1, f(x).diff(x).subs(x, x2): x3}`
*Details**
``eq`` can be any supported ordinary differential equation (see the
:py:mod:`~sympy.solvers.ode` docstring for supported methods).
This can either be an :py:class:`~sympy.core.relational.Equality`,
or an expression, which is assumed to be equal to ``0``.
``f(x)`` is a function of one variable whose derivatives in that
variable make up the ordinary differential equation ``eq``. In
many cases it is not necessary to provide this; it will be
autodetected (and an error raised if it couldn't be detected).
``hint`` is the solving method that you want dsolve to use. Use
``classify_ode(eq, f(x))`` to get all of the possible hints for an
ODE. The default hint, ``default``, will use whatever hint is
returned first by :py:meth:`~sympy.solvers.ode.classify_ode`. See
Hints below for more options that you can use for hint.
``simplify`` enables simplification by
:py:meth:`~sympy.solvers.ode.odesimp`. See its docstring for more
information. Turn this off, for example, to disable solving of
solutions for ``func`` or simplification of arbitrary constants.
It will still integrate with this hint. Note that the solution may
contain more arbitrary constants than the order of the ODE with
this option enabled.
``xi`` and ``eta`` are the infinitesimal functions of an ordinary
differential equation. They are the infinitesimals of the Lie group
of point transformations for which the differential equation is
invariant. The user can specify values for the infinitesimals. If
nothing is specified, ``xi`` and ``eta`` are calculated using
:py:meth:`~sympy.solvers.ode.infinitesimals` with the help of various
heuristics.
``ics`` is the set of boundary conditions for the differential equation.
It should be given in the form of ``{
f(x0): x1, f(x).diff(x).subs(x, x2):
x3}`` and so on. For now initial conditions are implemented only for
power series solutions of first-order differential equations which should
be given in the form of ``{
f(x0): x1}`` (See issue 4720). If nothing is
specified for this case ``f(0)`` is assumed to be ``C0`` and the power
series solution is calculated about 0.
``x0`` is the point about which the power series solution of a differential
equation is to be evaluated.
``n`` gives the exponent of the dependent variable up to which the power series
solution of a differential equation is to be evaluated.
例题2:求下列微分方程的解
from sympy.abc import x #引进符号变量x
from sympy import Function,diff,dsolve,sin
y=Function('y')
eq=diff(y(x),x,2)+2*diff(y(x),x)+2*y(x)-sin(x) #定义方程
con={
y(0): 0,diff(y(x),x).subs(x,0): 1}
y=dsolve(eq,ics=con)
print(simplify(y))
解出结果为:
E q ( y ( x ) , ( C 1 ∗ s i n ( x ) + C 2 ∗ c o s ( x ) ) ∗ e x p ( − x ) + s i n ( x ) / 5 − 2 ∗ c o s ( x ) / 5 ) Eq(y(x), (C1*sin(x) + C2*cos(x))*exp(-x) + sin(x)/5 - 2*cos(x)/5) Eq(y(x),(C1∗sin(x)+C2∗cos(x))∗exp(−x)+sin(x)/5−2∗cos(x)/5)
C1,C2代表任意常数
例题3: 求下列微分方程组的解:
method1:
import sympy as sp
t=sp.symbols('t')
x1,x2,x3=sp.symbols('x1,x2,x3',cls=sp.Function)
eq=[x1(t).diff(t)-2*x1(t)+3*x2(t)-3*x3(t),
x2(t).diff(t)-4*x1(t)+5*x2(t)-3*x3(t),
x3(t).diff(t)-4*x1(t)+4*x2(t)-2*x3(t)]
con={
x1(0):1, x2(0):2, x3(0):3}
s=sp.dsolve(eq,ics=con);
print(s)
method2(more simple one–matrix):
import sympy as sp
t=sp.symbols('t')
x1,x2,x3=sp.symbols('x1,x2,x3',cls=sp.Function)
x=sp.Matrix([x1(t),x2(t),x3(t)])
A=sp.Matrix([[2,-3,3],[4,-5,3],[4,-4,2]])
eq=x.diff(t)-A*x
s=sp.dsolve(eq,ics={
x1(0):1, x2(0):2, x3(0):3})
print(s)
The result:
[Eq(x1(t), C2*exp(-t) + C3*exp(2*t)), Eq(x2(t), C1*exp(-2*t) + C2*exp(-t) + C3*exp(2*t)), Eq(x3(t), C1*exp(-2*t) + C3*exp(2*t))]
c1,c2···代表任意常数,可以将初始值代回去得到一些常数具体的值