1.微分方程1.1基本概念微分方程是描述具有时间和空间的系统的状态演化的数学工具。物理学中许多涉及变力的运动学和动力学问题,如空气阻力是速度的函数的落体,都可以用微分方程求解。微分方程也广泛应用于化学、工程、经济学和人口统计学。具体来说,微分方程是指包含未知函数及其导数的关系表达式。
微分方程分为只有一个自变量的常微分方程和有两个或两个以上自变量的偏微分方程。微分方程按阶次可分为一阶、二阶和高阶。微分方程的阶取决于方程中最高导数的阶。微分方程还可以分为:(非)齐次、常(变)系数、(非线性)、初值问题/边界问题.1.2微分方程的数学建模微分方程的数学建模其实并不复杂。基本流程是分析问题属于哪一类问题,可以选择什么微分方程模型,然后如何利用现有的微分方程模型进行建模。在数学、力学、物理、化学等学科的课程中,会针对本学科的各种问题建立合适的数学模型。在中学课程中,各门学科的数学模型主要是线性或非线性方程,而在大学物理和各种专业课程中,越来越多的以微分方程描述的数学模型出现。数学建模中的微分方程问题通常是这些专业课中比较简单的模型。专业课教材上介绍一个模型的时候,往往解释得很详细。只要问题类型明确,数学模型选得好,建模和求解并不难。除此之外,还有很多现成的关于问题背景、适用范围、假设和求解过程的内容,可以在写论文的时候复制和参考。小白害怕的原因是当他看到微分方程时他感到害怕。第二,他缺乏专业背景。他不知道去哪里查资料,不会判断问题类型,不知道选择什么模型,不善于从题目内容中获取模型参数,不知道如何编程求解。于是,老师说,一看这就是XXX的问题,所以很明显,可以用XXX的模型。小白说,我们把问题b改一下,本系列将从一个简单的微分方程模型入手,重点介绍微分方程数值解的编程实现,通过分析问题,建立模型,帮助小白建立信心和动力。希望学完这个系列,你会发现微分方程模型是数学建模中最容易的题型:模型找教材,建模找例题,求解有套路,讨论有套路,论文成绩好。
1.3微分方程的数值解法在学习专业课的时候,经常会推导和求解微分方程的解析解。小白对微分方程模型的恐惧是从高等数学的“微分方程”开始形成的,并通过专业课程不断强化。实际上,只有少数微分方程可以解析求解,大多数微分方程只能用数值方法求解。微分方程的数值解法是将时间和空间离散化,然后微分成差分,建立递推关系,然后迭代计算得到任意时间和空间值。如果你还是觉得头晕,我们可以简化一下。建模就是把专业课课本上的公式抄下来,求解就是把公式的参数输入Python函数。先说解决。求解常微分方程的基本方法有欧拉法、龙格-库塔法等。详见各种教材。写数学模型竞赛论文还是可以抄几段的。本文遵循“编程方案”的概念,不涉及这些算法的具体内容。只讨论零基础如何使用Python的工具包和库函数求解微分方程模型。我们的选择是Python工具包三剑客:Scipy,
Scipy是Python算法库和数学工具包,包括最优化、线性代数、积分、插值、特殊函数、傅立叶变换、信号与图像处理、常微分方程求解等模块。介绍了Scipy是Python语言的Matlab,所以大部分数学建模问题都可以用它来解决。Numpy提供了线性代数运算、傅立叶变换、随机数生成等实现和计算高维数组的功能,还提供了与C/C等语言的集成工具。Matplotlib是一个可视化工具包,可以方便地绘制各种数据可视化图表,如折线图、散点图、直方图、条形图、箱线图、饼图、三维图等等。对了,还有一个Python符号运算工具包SymPy,解析求解积分和微分方程,也就是说,给出的结果都是微分方程的解析解表达式。很牛逼,但是只能用解析解解微分方程。所以,知道就好。
2.SciPy求解常微分方程(组)2.1一阶常微分方程(组)模型给定初始条件的一阶常微分方程(组)的标准型为:
公式中的y在常微分方程中是标量,在常微分方程中是数组向量。
2.2 SciPy.integrate.odeint()函数SciPy提供了两种求解常微分方程的方法:基于odeint函数的API相对简单易学,基于ode类的面向对象API更加灵活。**scipy.integrate.odeint() **是求解微分方程的具体方法,常微分方程采用数值积分求解。一阶刚体系统和非刚体系统的初值问题,可以用odeint函数中FORTRAN库odepack中的lsoda来解决。
scipy.integrate.odeint(func,y0,t,args=(),Dfun=None,col_deriv=0,full_output=0,ml=None,mu=None,rtol=None,atol=None,tcrit=None,h0=0.0,Hmax=0.0,hmin=0,ixpr=0,mxstep=0,mxhnil=0,mxordn=12,mxords=5,printmessg=0,tfirst=false)
Func:可调用(y,t,…)导数函数f (y,t),即y在t处的导数,表示为函数y0: array:初始条件y0。对于常微分方程y ^ 0,是数组向量t: array:解函数值对应的时间点的序列。序列的第一个元素是对应于初始条件y 0时间T0的初始值;时间序列必须是单调递增或单调递减的,并且允许重复值。其他参数简要介绍如下:
Args:将参数传递给func函数。当导数函数f (y,t,p 1,p 2,)包括可变参数P1、P2.参数P1,P2.可以通过args=(P1,P2,).有关argus的用法,请参见2.4中的示例2。DFUN:雅可比矩阵Dfun: func,行首。如果没有给定Dfun,则自动推导算法。Col_deriv:自动导出Dfun的方法。Printmessg:布尔值。控制是否打印收敛信息。其他参数用于控制求解算法的参数,一般可以忽略。int: y: array数组的主返回值,形状为(len(t),len(y0),给出时间的序列T中每一时刻的Y值。
3.例1: Scipy求解一阶常微分方程3.1例1:求微分方程的数值解3.2常微分方程的编程步骤
本文以此题为例,讲解scipy.integrate.odeint()求解常微分方程初值问题的步骤:1。导入scipy、numpy、matplotlib包;2.定义(t 2)中的导函数f (y,t)=s;3.定义初始值y 0和y的定义区间[t 0,t];4.调用odeint()求y在定义区间[t 0,t]内的数值解。
3.3 Python例程# 1。求解微分方程初值问题(scipy.integrate.odeint)从scipy.integrate导入odeint #导入scipy.integrate模块导入numpy as NP导入matplotlib . py plot as PLT def dy _ dt(y,t): #定义函数f(y,T)返回NP。sin (t * * 2) y0=[1] # y0=1或t=NP。arange (-10,10,0.01) # (start,stop,step) y=odeint (dy _ dt,y0,t) #解微分方程初值问题。y)工厂。标题(' scipy.integrate.odeint') PLT。Show () 3.4 Python例程运行结果4。例2: Scipy求解一阶常微分方程4.1例2:求Lorenz方程的数值解Lorenz混沌吸引子的轨迹可以用下面三个微分方程来描述:
洛伦兹方程将大气流体运动的强度X与水平和垂直温度变化Y和Z联系起来,模拟大气对流系统,已广泛应用于天气预报、空气污染和全球气候变化的研究中。参数\a称为plante数,\为归一化瑞利数,\与几何形状有关。洛伦兹方程是一组非线性微分方程,不能解析求解,只能用数值方法求解。
4.2洛伦兹方程问题的编程步骤以此题为例讲解scipy.integrate.odeint()解决常微分方程初值问题的步骤:1。导入scipy、numpy、matplotlib包;2.定义导数函数lorenz(W,t,p,r,b)。注意odeint()函数中定义导函数的标准形式是f (y,t),对于微分方程,y代表向量。为避免混淆,我们记为W=[x,y,z],函数lorenz(W,t)定义了导函数f (W,t)。p、R、B、R和B用于表示方程中的参数、和,导数定义函数编程如下:
#导数函数,求W=[x,Y,z]DX _ DT=P *(Y-X)# DX/DT P:sigma dy _ DT=X *(R-z)-Y # dy/DT=X *(R-z)-Y,R:rho dz _ DT=X * Y-B * z # dz/DT=X * Y-B * z,B;Beta返回NP。数组([DX _ DT,DY _ DT,DZ _ DT]) 3。定义初始值W 0和W的定义区间[T0,T];4.调用odeint()求W WW在定义区间[t 0,t]内的数值解。注意,在例程中,参数(p,r,b)通过args=paras或args=(10.0,28.0,3.0)传递给导数函数lorenz(W,t,p,r,b)。当然,参数(p,r,b)可以直接在导函数lorenz()中设置,而不是作为函数参数传递。而例程的参数传递法使导函数的结构更清晰,更通用。另外,对于变参数问题,使用这种参数传递方法非常方便。
4.3洛伦兹方程问题# 2的Python例程。解决微分方程初值问题。integrate.odeint)来自scipy。integrate Import ode int # Import scipy . integrate module Import numpy as NP Import matplotlib . py plot as PLT from mpl _ toolkits . mplot 3d Import axes 3d #导数函数,求W=[x,Y,导数DW/DT DEF Lorenz (W,T,P,R,B) of point z]: # byYoucansX,Y,Z=W # W=[X,Y,Z] DX _ DT=P * (Y-X) # DX/DT=P * (Y-X) pBETA返回NP。数组([DX _ DT,DY _ DT,DZ _ DT]) T=NP。ARANGE (0,30,0.01)# Created时间 points(start,stop,step) PARAS=(10.0,28.0,3.0) R,b)#调用ode求解lorenz,并使用两个不同的初始值W1,W2求解W1=(0.0,1.00,0.0) #将初始值定义为W1Track1=odeint (Lorenz,W1,T,Args=(10.0,28.0,3。 Args=paras) #通过para传递导数函数的参数# Drawing Figure=PLT . Figure()ax=axes3d(Figure)ax . plot(track1 [:0],track1 [:1],track 1[:2],Color=' magenta ')# Drawing track 1ax . plot(track2 [:0],track2 [:1],track 2[:2],Color=' deep sky blue ')# Drawing track 2ax . set _ title(' Lorenz Attractor by scipy . integrate . odeint '
4.4洛伦兹方程问题Python例程运行结果5。例3: Scipy必须替换变量求解高阶常微分方程,变成一阶微分方程,然后用odeint求数值解。
5.1例3:求解二阶RLC振荡电路的数值解输入响应为零的RLC振荡电路可以用下面的二阶微分方程来描述:
制造
当零输入响应u s=0时,上述公式可以写成:
对于二阶微分方程问题,引入变量v=d u/d t v,通过变量替换将原方程转化为如下微分方程组:
这样我们就可以用上一节解微分方程的方法来解高阶微分方程了。
5.2二阶微分方程的编程步骤:以RLC振荡器电路为例讲解scipy.integrate.odeint()求解高阶常微分方程初值问题的步骤:1 .导入scipy、numpy、matplotlib包;2.定义导函数deriv(Y,t,a,w)。注意odeint()函数中定义导函数的标准形式是f (y,t)。在这个问题中,Y代表一个向量,记为Y=[u,v]。导数定义函数deriv(Y,t,a,w)编程如下,其中a和w分别表示方程中的参数。
#导数函数,求该点的导数dy/dt def deriv (y,t,a,w)Y=[u,v] dy _ dt=[v,-2 * a * v-w * w * u]返回dy _ dt3。4.调用odeint()求Y=[u,v]在定义区间[t 0,t]内的数值解。在例程中,通过args=paras将参数(a,w)传递给导数函数deriv(Y,t,a,w)。在本例中,我们将检查不同参数对结果的影响。这种参数传递方法使用起来非常方便。
5.3二阶微分方程问题的Python例程# 3。求解二阶微分方程初值问题(scipy.integrate . odeint)# Secondode by scipy . integrate . odeint from scipy . integrate Import odeint #导入scipy . integrate模块的导数函数导入numpy as NP导入matplotlib.py plot as PLT #,求导数dy/dt def deriv (y,t,a,w): u,v=y # y=[u,V] dy _ dt=[v,-2 * a * V-w * w * w * u]返回dy _ dtt=np0.6) #过阻尼:a^2-w^2 0paras2=(1,1) #临界阻尼:a^2-w^2=0paras3=(0.3,1) #欠阻尼:一个2-W20 #调用ode解决问题,用两个不同的初始值W1,0.0) #定义初始值为y0=[u0,v0] y1=odeint (deriv,y0,T,args=paras1) # args设置参数Y2=odeint(deriv,y0,T,Args=paras2) # args设置参数y3=odeint的导数函数(deriv,Y0,T,args=paras3) # args设置参数# W2=(0.0,1.01,0.0) #定义初始值为W2# track2=odeint args=paras) #通过para #传递导数函数的参数绘图PLT.plot (t,y1 [:0],' r-'label=' u1 (t)') PLT.plot (t,y2 [:0],' b-'label='' 'g-'label='u3(t)')plt.plot(t,Y1[:1],' r:'label='v1(t)')plt.plot(t,Y2[:1],' b:'
5.4二阶方程问题Python例程运行结果讨论:RLC串联电路是典型的二阶系统。在零输入条件下,根据\和\的关系,电路的输出响应有四种情况:
过阻尼: 2 20,有两个不相等的负实根;临界阻尼: 2 2=0,有两个相等的负实根;欠阻尼: 2 20,有一对共轭复根;无阻尼:R=0 R=0R=0,有一对纯虚根。程序中选取的三组参数分别对应过阻尼、临界阻尼和欠阻尼的情况,微分方程的数值结果很好地反映了不同情况对应的曲线。
6.总结小白首先要有信心。利用Scipy工具包求解标准形式的微分方程(组),编程实现非常简单,容易掌握。其次要认识到,由于微分方程解的特征与模型参数密切相关,不同参数的解可能具有完全不同的形状。这就涉及到对模型的稳定性和敏感性的分析,很容易把论文写得非常丰富精彩。不能从问题出发建立微分方程模型,或者讨论参数对稳定性和灵敏度的影响,怎么办?谁让你自己做的?当然,你要先找相关专业的教材和论文,从中选择相对接近、简单的理论和模型,然后通过各种假设,把题目强行简化到模型中的条件,这样就可以画虎为猫了。这一部分的结束。】58660 . 68686868661
版权声明:本文为CSDN博主“youcans”原创文章,遵循CC 4.0 BY-SA版权协议。请附上原始文章和本声明来源的链接。原文链接:https://blog.csdn.net/youcans/article/details/117702996