Python求解常微分方程的符号解+数值解例题

少年开新

🌾常微分方程的符号解Python提供了功能强大的求解常微分方程符号解的函数dsolve。下面我们来看几道例题。🌼例1求解微分方程y′=−2y+2x2+2x,y(0)=1y^{'}=-2y+2x^2+2x,y(0)=1y′=−2y+2x2+2x,y(0)=1代码如下：importsympyasspx=sp.var('x')y=sp.Function('y')eq=y(x).diff(x)+2*y(x)-2*x*x-2*xs=sp.dsolve(eq,ics={y(0):1})s=sp.simplify(s)print(s)1234567运行结果如下：在交互式窗口中运行得到：🌼例2求解二阶微分方程y′′−2y′+y=ex,y(0)=1,y′(0)=−1y^{''}-2y^{'}+y=e^x,y(0)=1,y'(0)=-1y′′−2y′+y=ex,y(0)=1,y′(0)=−1代码如下：importsympyasspx=sp.var('x')y=sp.Function('y')eq=y(x).diff(x,2)-2*y(x).diff(x)+y(x)-sp.exp(x)con={y(0):1,y(x).diff(x).subs(x,0):-1}s=sp.dsolve(eq,ics=con)print(s)1234567解得方程为：🌼例3已知输入信号为u(t)=e−tcos(t)u(t)=e^{-t}cos(t)u(t)=e−tcos(t),试求下面微分方程的解。y(4)(t)+10y(3)(t)+35y′′(t)+50y′(t)+24y(t)=u′′(t)y^{(4)}(t)+10y^{(3)}(t)+35y^{''}(t)+50y^{'}(t)+24y(t)=u^{''}(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)=1y(0)=0,y^{'}(0)=-1,y^{''}(0)=1,y^{'''}(0)=1y(0)=0,y′(0)=−1,y′′(0)=1,y′′′(0)=1代码如下：importsympyasspt=sp.var('t')y=sp.Function('y')u=sp.exp(-t)*sp.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=sp.dsolve(eq,ics=con)s=sp.expand(s)print(s)123456789求得解如下：🌾常微分方程的数值解🌼例4求解例1微分方程y′=−2y+2x2+2x,y(0)=1y^{'}=-2y+2x^2+2x,y(0)=1y′=−2y+2x2+2x,y(0)=1的数值解。代码：fromscipy.integrateimportodeintimportnumpyasnpimportpylabaspltimportsympyasspdy=lambday,x:-2*y+2*x*x+2*xxx=np.linspace(0,3,31)#自变量的取值s=odeint(dy,1,xx)#y的取值print('x={}\n对应的数值解y={}'.format(xx,s.flatten()))plt.plot(xx,s,'*')x=sp.var('x')y=sp.Function('y')eq=y(x).diff(x)+2*y(x)-2*x*x-2*xs2=sp.dsolve(eq,ics={y(0):1})sx=sp.lambdify(x,s2.args[1],'numpy')plt.plot(xx,sx(xx))#观察吻合度plt.show()12345678910111213141516171819运行结果：🌼例5求解微分方程(1−x)d2ydx2=151+(dydx)2(1-x)\frac{d^2y}{dx^2}=\frac{1}{5}\sqrt{1+(\frac{dy}{dx})^2}(1−x)dx2d2y​=51​1+(dxdy​)2​0＜x≤1的数值解。求数值解时，需要把二阶微分方程转化为一阶微分方程，引进y1=y,y2=y′y_1=y,y_2=y^{'}y1​=y,y2​=y′,则可以将方程转化为如下的一阶微分方程组：y1′=y2,y1(0)=0y_1^{'}=y_2,y_1(0)=0y1′​=y2​,y1​(0)=0y2′=15(1−x)1+y22,y2(0)=0y_2^{'}=\frac{1}{5(1-x)}\sqrt{1+y_2^2},y_2(0)=0y2′​=5(1−x)1​1+y22​​,y2​(0)=0代码：fromscipy.integrateimportodeintimportnumpyasnpimportpylabaspltyx=lambday,x:[y[1],np.sqrt(1+y[1]**2)/5/(1-x)]x0=np.arange(0,1,0.00001)y0=odeint(yx,[0,0],x0)plt.rc('font',size=16)plt.plot(x0,y0[:,0])plt.show()12345678910运行结果：