数模竞赛-常微分方程与差分方程

常见模型

略。

常微分方程求解-符号解

例题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} #subs替换x为0
s=sympy.dsolve(eq,ics=con)
print(s)

例题3：

已知输入信号为$u(t)=e^{-t}\cos 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：

求：

其中：

代码：

1
2
3
4
5
6
7
8
import sympy
sympy.var('t')
sympy.var('x1:4',cls=sympy.Function) #定义3个符号函数
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) #Runge-Kutta法
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：求下式的数值解：

先引进$y_1=y,y_2=y^{‘}$，化为一阶微分方程：

代码：

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]) #0表示求第一个方程的解
pylab.show()

例题3：

洛伦兹模型的混沌效应：

其中$\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()