matlab-常微分方程数值解法 1

科学计算与MATLAB，讲座：2012，唐建国中南大学材料科学与工程学院，第7讲常微分方程数值解，小结，介绍欧拉逼近法龙格-库塔(R-K)法，MATLAB常微分方程函数小结，1。引言，物理学研究的许多物质运动过程都是用常微分方程来描述的。例如，粒子的加速运动、简谐振动等。简单的问题可以用解析方法解决，而大多数实际问题则用数值方法解决。一阶常微分方程初值问题：数值解法是在一些离散节点xn上求y(x)的近似值yn来逼近Y (xn)，二阶常微分方程，2.1简单欧拉(L.Euler，1707-1783)方法。欧拉数值算法是通过递归求解初始值。递归解从初始值开始，后一个函数值从前一个函数值获得。关键是构造递推公式。2，欧拉近似法，欧拉数值算法递归公式构造，2.1.1差分法，差分法是用差商近似代替导数，即用微分方程代替求积，对于等距节点，可以求积，在xn节点上，微分方程可以写成，作如下近似，得到欧拉解递归公式的一般形式：具体求解过程是：简单欧拉方法程序函数outx，OUTY=MYUELER(FUN，x0，XT，Y0， PointNum)% MYEULER使用前向差分EULER方法求解微分方程%fun表示f(x，y)%x0，xt表示自变量的初始值和终值%y0表示函数在x0处的值，可以用向量形式%PointNum表示自变量在x0，xt上取的点数5|PointNum=0%。 如果该函数仅输入4个参数值，则POINTNUM默认为100PointNum=100endifnargindsolve(dy 3*x*y=x*exp(-x2)ans=(1/3*exp(-x*(x-3*t)C1)* exp(-3 * x * t)dsolve(dy 3*x*y=x*exp(-x2),x)ans=exp(-x2)exp(-3/2*x2)*c1，示例：使用MATLAB的符号解来求解常微分方程的特殊解：dSolve (x * dy2 * y-exp (x)=0，y (1)=2 * exp (1)，x)ans=(exp(x)* x-exp(x)2 * exp(1)/x 2，示例：ODE45用于求解描述刚性问题的微分方程：functiondy=刚性(1)%行向量dy(1)=y(2)* y(3)；dy(2)=-y(1)* y(3)；dy(3)=-0.51 * y(1)* y(2)；options=odeset(RelTol，1e-4，AbsTol，1e-41e-41e-5)；T，Y=ode45(刚性，012，011，选项)；图(t，y (:1)，-，t，y (:2)，-。t，y (:3)，)图例(y1，y2，y3)，示例：使用MATLAB函数ode23，ode45，ode113求解多阶常微分方程：functiondy=myfun03 (x，y)dy=0(3，1)%初始化变量dydy(1)=y(2)；%dy(1)表示y的一阶导数，它等于y的第二列值dy(2)=y(3 );%dy(2)表示二阶导数dy(3)=2 * y(3)/x 3 3 * y(2)/x 3 3 * exp(2 * x)/x 3% dy(3)表示y的三阶导数。x45，y45=ode45(myfun03，1，10，11030)；x113，y113=ode113(myfun03，1，10，11030)；图(1)%第一个图形的绘图(X23，Y23 (:1)，* R，x45，Y45 (:1)，OB，X113，Y113 (:1)，G)%制作图例(ode23解决方案，ode45解决方案，Ode113解决方案)图(2)绘图(X45，y45)%以ode45为例制作函数及其导数图形图例(Y，Y一阶导数，Y二阶导数)标题(Y，Y二阶导数)计算过程中涉及许多计算方法。虽然热传导是一种易于数学处理的传热方法，但其过程往往涉及常微分、偏微分方程和线性(非线性)方程的求解。有一个内部冷却的钢导体，其外径为4厘米，内径为1.5厘米，电流密度为5000安培/平方厘米。导体单位时间发出的热量等于流体同时带走的热量，导体内壁表面温度为70。假设外壁表面完全绝缘。尝试确定导体内部的温度分布。已知钢的热导率k=0.38Kw/mK，电导率=21011km。解决方法：这是一维的在命令窗口中执行以下代码：cleard _ equi=D2t Dt/r 131579=0。条件=(t0.0075)=70，D(t 0.02)=0 ；%边界条件t=d求解(d _ equal，condition，r)程序执行结果：程序执行结果：即温度分布方程为：工程中遇到的热传导问题，由于物体的几何形状复杂或边界条件难以描述，往往无法用数值方法解决。常微分方程包括两类初值问题和边值问题。常通过调用ode45()和ode23()函数来实现初值问题的数值解。在MATLAB编辑器中编写函数BZ.m并保存。functionBZclearallclca=0.0075b=0.02%r值的范围solinit=BVP ini(t line space(a，b，10)，7080)。%使用初始值对解决方案初始化SOL=BVP4C ( ODEFEN， BC FUN，SOLIIT)；%求解ODE方程r=0.007533600.0012533600.02t=deva(LSOL，r)函数dtdr=ODEFN (r，t)%定义ODE方程函数dtdr=(T2)；-131580-1/r *(T2)；Function