数模竞赛-常微分方程与差分方程
常见模型
略。
常微分方程求解-符号解
例题1：
求解$y^{‘}=-2y+2x^2+2x,y(0)=1$：
1
2
3
4
5
6
7
| import sympy
sympy.var('x')
y=sympy.Function('y')
eq=y(x).diff(x)+2*y(x)-2*x**2-2*x
s=sympy.dsolve(eq,ics={y(0):1})
s=sympy.simplify(s)
print(s)
|
得$y=x^2+e^{-2x}$。
例题2：
求解$y^{‘’}-2y^{‘}+y=e^x,y(0)=1,y^{‘}(0)=-1$：
1
2
3
4
5
6
7
| import sympy
sympy.var('x')
y=sympy.Function('y')
eq=y(x).diff(x,2)-2*y(x).diff(x)+y(x)-sympy.exp(x)
con={y(0):1,y(x).diff(x).subs(x,0):-1}
s=sympy.dsolve(eq,ics=con)
print(s)
|
例题3：
已知输入信号为$u(t)=e^{-t}\cos t$，求解：
$$
y^{(4)}(t)+10y^{(3)}(t)+35y^{‘’}(t)+50y^{‘}(t)+24y(t)=u^{‘’}(t)
$$
其中$y(0)=0,y^{‘}(0)=-1,y^{‘’}(0)=1,y^{‘’’}(0)=1$。
代码：
1
2
3
4
5
6
7
8
9
| import sympy
sympy.var('t')
y=sympy.Function('y')
u=sympy.exp(-t)*sympy.cos(t)
eq=y(t).diff(t,4)+10*y(t).diff(t,3)+35*y(t).diff(t,2)+50*y(t).diff(t)+24*y(t)-u.diff(t,2)
con={y(0):0,y(t).diff(t).subs(t,0):-1,y(t).diff(t,2).subs(t,0):1,y(t).diff(t,3).subs(t,0):1}
s=sympy.dsolve(eq,ics=con)
s=sympy.expand(s)
print(s,s.args[1])
|
例题4：
求：
$$
\begin{cases}
\displaylines{
\dfrac{\mathrm{d}x}{\mathrm{d}t}=Ax,\\
x(0)=[1,1,1]^{\mathrm{T}}
}
\end{cases}
$$
其中：
$$
x(t)=\left[x_1(t),x_2(t),x_3(t)\right]^{\mathrm{T}},A=\left[
\begin{matrix}
\displaylines{
3&-1&1\\
2&0&-1\\
1&-1&2
}
\end{matrix}
\right]
$$
代码：
1
2
3
4
5
6
7
8
| import sympy
sympy.var('t')
sympy.var('x1:4',cls=sympy.Function)
x=sympy.Matrix([x1(t),x2(t),x3(t)])
A=sympy.Matrix([[3,-1,1],[2,0,-1],[1,-1,2]])
eq=x.diff(t)-A@x
s=sympy.dsolve(eq,ics={x1(0):1,x2(0):1,x3(0):1})
print(s)
|
常微分方程求解-数值解
数值解即为符号解的自变量代入具体数值。
例题1：求$y^{‘}=-2y+2x^2+2x,y(0)=1$的符号解和数值解。
1
2
3
4
5
6
7
8
9
10
11
12
13
| import scipy,numpy,pylab,sympy
dy=lambda y,x:-2*y+2*x**2+2*x
xx=numpy.linspace(0,3,31)
s=scipy.integrate.odeint(dy,1,xx)
print(xx,s.flatten())
pylab.plot(xx,s,'*')
sympy.var('x')
y=sympy.Function('y')
eq=y(x).diff(x)+2*y(x)-2*x**2-2*x
s2=sympy.dsolve(eq,ics={y(0):1})
sx=sympy.lambdify(x,s2.args[1],'numpy')
pylab.plot(xx,sx(xx))
pylab.show()
|
例题2：求下式的数值解：
$$
\begin{cases}
\displaylines{
(1-x)\dfrac{\mathrm{d}^2y}{\mathrm{d}x^2}=\dfrac{1}{5}\sqrt{1+\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right)^2},&0<x\leqslant1,\\
y(0)=0,y^{‘}(0)=0.
}
\end{cases}
$$
先引进$y_1=y,y_2=y^{‘}$，化为一阶微分方程：
$$
\begin{cases}
\displaylines{
y_1^{‘}=y_2,&y_1(0)=0,\\
y_2^{‘}=\dfrac{1}{5(1-x)}\sqrt{1+y_2^2},&y_2(0)=0.
}
\end{cases}
$$
代码：
1
2
3
4
5
6
7
| import scipy,numpy,pylab
yx=lambda y,x:[y[1],numpy.sqrt(1+y[1]**2)/5/(1-x)]
x0=numpy.arange(0,1,0.00001)
y0=scipy.integrate.odeint(yx,[0,0],x0)
pylab.rc('font',size=16)
pylab.plot(x0,y0[:,0])
pylab.show()
|
例题3：
洛伦兹模型的混沌效应：
$$
\displaylines{
\begin{cases}
\dot{x}=\sigma(y-x),\
\dot{y}=\rho x-y-xz,\
\dot{z}=xy-\beta z
\end{cases}
}
$$
其中$\sigma$为普朗特数，$\rho$为瑞利数，$\beta$为方向比。
代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
| import scipy,numpy,pylab
numpy.random.seed(2)
sigma=10
rho=28
beta=8/3
g=lambda f,t:[sigma*(f[1]-f[0]),rho*f[0]-f[1]-f[0]*f[2],f[0]*f[1]-beta*f[2]]
s01=numpy.random.rand(3)
t0=numpy.linspace(0,50,5000)
s1=scipy.integrate.odeint(g,s01,t0)
pylab.rc('text',usetex=True)
pylab.rc('font',size=16)
pylab.subplots_adjust(wspace=0.6)
ax=pylab.subplot(121,projection='3d')
pylab.plot(s1[:,0],s1[:,1],s1[:,2],'r')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
s02=s01+0.000001
s2=scipy.integrate.odeint(g,s02,t0)
pylab.subplot(122)
pylab.plot(t0,s1[:,0]-s2[:,0],'.-')
pylab.xlabel('$t$')
pylab.ylabel('$x_1(t)-x_2(t)$',rotation=90)
pylab.show()
|
