Matlab基础：第6讲 求微分方程的解

《Matlab基础：第6讲 求微分方程的解》由会员分享，可在线阅读，更多相关《Matlab基础：第6讲 求微分方程的解（30页珍藏版）》请在装配图网上搜索。

1、第六讲求微分方程的解Matlab基础q 自牛顿发明微积分以来，微分方程在描述事物运自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。数学建模所得到的方程，绝大多数是微分方程。q 由于实际应用的需要，人们必须求解微分方程。由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。数微分方程需要利用数值方法来近似求解。q 本节主要介绍求解微分方程（组）的本节主要介绍求解

2、微分方程（组）的 Matlab 命令命令.研究如何用研究如何用 Matlab 来计算微分方程（组）的数值解，来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的基本数值解法并重点介绍一个求解微分方程的基本数值解法Euler折线法折线法。问题背景和课程要求问题背景和课程要求求解微分方程求解微分方程(组组)的的Matlab 命令命令q 用用 Maltab自带函数自带函数 求解求解u 求解析解：求解析解：dsolve 可求通解可求通解,也可求特解也可求特解u 求数值解：求数值解：ode45、ode23、ode113、ode23t、ode15s、ode23s、ode23tbdsolve 求解析解

3、求解析解q dsolve 的使用的使用y=dsolve(eq1,eq2,.,cond1,cond2,.,v)其中其中 y 为输出，为输出，eq1、eq2、.为微分方程，为微分方程，cond1、cond2、.为初值条件，为初值条件，v 为自变量。为自变量。例例 1：求微分方程求微分方程 的通解，并验证。的通解，并验证。22xdyxyxedx y=dsolve(Dy+2*x*y=x*exp(-x2),x)syms x;diff(y)+2*x*y-x*exp(-x2)dsolve 的使用的使用q 几点说明几点说明l 如果省略初值条件，则表示求通解；如果省略初值条件，则表示求通解；l 如果省略自变量，

4、则默认自变量为如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x);dy/dx=2xdsolve(Dy=2*x);dy/dt=2xl 若找不到解析解，则返回其积分形式。若找不到解析解，则返回其积分形式。l 微分方程中用微分方程中用 D 表示对表示对 自变量自变量 的导数，如：的导数，如：Dy y；D2y y；D3y ydsolve 举例举例例例 2：求微分方程求微分方程 在初值条件在初值条件 下的特解，并画出解函数的图形。下的特解，并画出解函数的图形。0 xxyye y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x)ezplot(y);12()y

5、e v dsolve 举例举例例：例：x,y=dsolve(Dx+5*x=0,Dy-3*y=0,.x(0)=1,y(0)=1,t)r=dsolve(Dx+5*x=0,Dy-3*y=0,.x(0)=1,y(0)=1,t)这里返回的这里返回的 r 是一个是一个 结构类型结构类型 的数据的数据r.x%查看解函数查看解函数 x(t)r.y%查看解函数查看解函数 y(t)只有很少一部分微分方程（组）能求出解析解。只有很少一部分微分方程（组）能求出解析解。大部分微分方程（组）只能利用大部分微分方程（组）只能利用数值方法数值方法求数值解。求数值解。dsolve的输出个数只能为一个的输出个数只能为一个 或或

6、与方程个数相等。与方程个数相等。dsolve 举例举例例例3：求微分方程组求微分方程组 在初值条件在初值条件 下的特解，并画出解函数的图形。下的特解，并画出解函数的图形。530tdxxyedtdyxydt x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0,.x(0)=1,y(0)=0,t)ezplot(x,y,0,1.3);0010|ttxy 注：解微分方程组时，如果所给的输出个数与方程个数相同，注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解则方程组的解按词典顺序按词典顺序输出；如果只给一个输出，则输出输出；如果只给一个输出，则输出的是一个包含解的

7、的是一个包含解的结构结构（structure）类型的数据。类型的数据。微分方程的数值解微分方程的数值解。的相应近似值的相应近似值求出准确值求出准确值，值处，即对值处，即对开始的若干离散的开始的若干离散的其数值解是指由初始点其数值解是指由初始点，：对常微分方程对常微分方程nnnyyxyxyxxxxxy,y)(,),(),y(x x )y(xy)f(x,y 2121210000 q 常微分方程数值解的定义常微分方程数值解的定义q 考虑一维经典考虑一维经典初值问题初值问题00,(,)(),dyf x yy xyxa bdx u 基本思想：基本思想：用差商代替微商用差商代替微商根据根据 Talyor

8、公式，公式，y(x)在点在点 xk 处有处有2()()()()()kkky xy xxxyxOx 1kkhxx 11()()()()()kkkkky xy xy xy xdyO hdxhhx 21()()()()kkky xy xhyxO h Euler 折线法折线法初值问题的初值问题的Euler折线法折线法q 具体步骤：具体步骤：等距剖分：等距剖分：0121nnaxxxxxb 步长：步长：1 0 1 21()/,kkhxxknnba u 分割求解区间分割求解区间u 差商代替微商差商代替微商1()()kkky xy xdydxhx 1()()()kkky xy xh yx 得方程组：得方程组：

9、0011 ()(,)kkkkkkyy xyyh f xyxxh 分割求解区间，差商代替微商，解代数方程分割求解区间，差商代替微商，解代数方程 为分割点为分割点 0nkkx k=0,1,2,.,n-1yk 是是 y(xk)的近似的近似Euler 折线法举例折线法举例例：例：用用 Euler 法解初值问题法解初值问题22 0 201 ,()dyxydxyxy 取步长取步长 h=(2-0)/n=2/n，得差分方程得差分方程00110 1 2,(,)(/)kkkkkkkkkkxyyyh f xyyh yxyxxh 当当 h=0.4，即即 n=5 时，时，Matlab 源程序见源程序见 fulu1.m解

10、：解：Euler 折线法源程序折线法源程序clearf=sym(y+2*x/y2);a=0;b=2;h=0.4;n=(b-a)/h+1;%n=(b-a)/h;x=0;y=1;szj=x,y;for i=1:n-1%i=1:n y=y+h*subs(f,x,y,x,y);x=x+h;szj=szj;x,y;endszjplot(szj(:,1),szj(:,2),or-)Euler折线法举例（续）折线法举例（续）解析解：解析解：1 3352233/xyex 解析解解析解近似解近似解Runge-Kutta 方法方法q 为了减小误差，可采用以下方法：为了减小误差，可采用以下方法：u 让步长让步长 h

11、 取得更小一些；取得更小一些；u 改用具有较高精度的数值方法：改用具有较高精度的数值方法：q 龙格龙格-库塔方法库塔方法Runge-Kutta(龙格龙格-库塔库塔)方法方法u 是是一类一类求解常微分方程的数值方法求解常微分方程的数值方法u 有多种不同的迭代格式有多种不同的迭代格式Runge-Kutta 方法方法q 用得较多的是用得较多的是 四阶四阶R-K方法方法（教材第（教材第 79 页）页）00111234 (22)/6(),kkkkyy xxxhyyh LLLL 12132432222(,)(/,/)(/,/)(,)kkkkkkkkLf xyLf xhyhLLf xhyhLLf xh yh

12、L 其中其中四阶四阶 R-K 方法方法源程序源程序clear;f=sym(y+2*x/y2);a=0;b=2;h=0.4;n=(b-a)/h+1;%n=(b-a)/h;x=0;y=1;szj=x,y;for i=1:n-1%i=1:n l1=subs(f,x,y,x,y);l2=subs(f,x,y,x+h/2,y+l1*h/2);l3=subs(f,x,y,x+h/2,y+l2*h/2);l4=subs(f,x,y,x+h,y+l3*h);y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=szj;x,y;endplot(szj(:,1),szj(:,2),dg-)Run

13、ge-Kutta 方法方法Euler 法与法与 R-K法误差比较法误差比较Matlab函数数值求解函数数值求解T,Y=solver(odefun,tspan,y0)其中其中 y0 为初值条件，为初值条件，tspan为求解区间；为求解区间；Matlab在数值求解在数值求解时时自动对求解区间进行分割自动对求解区间进行分割，T(向量向量)中返回的是分割点的中返回的是分割点的值值(自变量自变量)，Y(向量向量)中返回的是解函数在这些分割点上的中返回的是解函数在这些分割点上的函数值。函数值。solver 为为Matlab的的ODE求解器求解器（可以是（可以是 ode45、ode23、ode113、ode

14、15s、ode23s、ode23t、ode23tb）没有一种算法可以有效地解决所有的没有一种算法可以有效地解决所有的 ODE 问题，因此问题，因此MATLAB 提供了多种提供了多种ODE求解器求解器，对于不同的对于不同的ODE，可以调用不同的可以调用不同的求解器求解器。Matlab提供的提供的ODE求解器求解器求解器求解器 ODE类型类型特点特点说明说明ode45非刚性非刚性单步法；单步法；4，5 阶阶 R-K 方法；方法；累计截断误差为累计截断误差为(x)3大部分场合的大部分场合的首选方法首选方法ode23非刚性非刚性单步法；单步法；2，3 阶阶 R-K 方法；方法；累计截断误差为累计截断误

15、差为(x)3使用于精度较低的情形使用于精度较低的情形ode113非刚性非刚性多步法；多步法；Adams算法；高低精算法；高低精度均可到度均可到 10-310-6计算时间比计算时间比 ode45 短短ode23t适度刚性适度刚性 采用梯形算法采用梯形算法适度刚性情形适度刚性情形ode15s刚性刚性多步法；多步法；Gears 反向数值微反向数值微分；精度中等分；精度中等若若 ode45 失效时，可失效时，可尝试使用尝试使用ode23s刚性刚性单步法；单步法；2 阶阶Rosebrock 算算法；低精度法；低精度当精度较低时，计算时当精度较低时，计算时间比间比 ode15s 短短ode23tb刚性刚性

16、梯形算法；低精度梯形算法；低精度当精度较低时，计算时当精度较低时，计算时间比间比ode15s短短参数说明参数说明odefun 为为显式常微分方程显式常微分方程，可以用命令，可以用命令 inline 定义，或定义，或在在函数文件函数文件中定义，然后通过函数句柄调用。中定义，然后通过函数句柄调用。fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);注：注：也可以在也可以在 tspan 中指定对求解区间的分割，如：中指定对求解区间的分割，如：x,y=ode23(fun,0:0.1:0.5,1);%此时此时 x=0:0.1:0.5T,Y=solve

17、r(odefun,tspan,y0)求初值问题求初值问题 的数值解，求解范的数值解，求解范围为围为 0,0.5222201()dyyxxdxy 例例 4：数值求解举例数值求解举例如果需求解的问题是如果需求解的问题是高阶高阶常微分方程，则需将其化为常微分方程，则需将其化为一阶常一阶常微分方程组微分方程组，此时需用，此时需用函数文件函数文件来定义该常微分方程组。来定义该常微分方程组。122212112101 00 7/()(),(),dxdtxdxdtxxxxx 令令 ，则原方程可化为，则原方程可化为12,dyxy xdt 求解求解 Ver der Pol 初值问题初值问题2221001 00 7

18、()(),(),d ydyyydtdtyy 例例 5：数值求解举例数值求解举例l 先编写函数文件先编写函数文件 verderpol.mfunction xprime=verderpol(t,x)global mu;xprime=x(2);mu*(1-x(1)2)*x(2)-x(1);l 再编写脚本文件再编写脚本文件 vdpl.m，在命令窗口直接运行该文件。，在命令窗口直接运行该文件。clear;global mu;mu=7;y0=1;0;t,x=ode45(verderpol,0,40,y0);plot(t,x(:,1),r-);Matlab 求解微分方程小结求解微分方程小结q Matlab

19、函数函数u 求解析解（求解析解（通解通解或特解），用或特解），用 dsolveu 求数值解（特解），用求数值解（特解），用 ode45、ode23.q Matlab 编程编程u Euler 折线法折线法u Runga-Kutta 方法方法上机作业上机作业u 求微分方程求微分方程xeyyyxsin52 的通解的通解00yxdtdyyxdtdx0|,1|00ttyx在初始条件在初始条件下的特解，并画出解函数的图形下的特解，并画出解函数的图形u 求微分方程组求微分方程组1)0(,cosyxeyyx的数值解(步长h取0.1)，求解范围为区间0,3u用用Euler法及四阶法及四阶 Runge-Kutta

20、 法求解微分方程初值法求解微分方程初值问题问题作业题解答作业题解答l 第一题第一题 cleary=dsolve(D2y-2*Dy+5*y=exp(x)*sin(x),x)y=exp(x)*sin(2*x)*C2+exp(x)*cos(2*x)*C1+1/3*exp(x)*sin(x)作业题解答作业题解答l 第二题第二题 clearx,y=dsolve(Dx+x+y=0,Dy+x-y=0,x(0)=1,y(0)=0,t)x=(-1/4*2(1/2)+1/2)*exp(2(1/2)*t)+(1/4*2(1/2)+1/2)*exp(-2(1/2)*t)y=-(-1/4*2(1/2)+1/2)*2(1

21、/2)*exp(2(1/2)*t)+(1/4*2(1/2)+1/2)*2(1/2)*exp(-2(1/2)*t)-(-1/4*2(1/2)+1/2)*exp(2(1/2)*t)-(1/4*2(1/2)+1/2)*exp(-2(1/2)*t)作业题解答作业题解答l 第二题第二题 ezplot(x,y,0,1.0);00.20.40.60.811.21.41.6-1.2-1-0.8-0.6-0.4-0.20 xyx=(-1/4 21/2+1/2)exp(21/2 t)+(1/4 21/2+1/2)exp(-21/2 t),y=-.-(1/4 21/2+1/2)exp(-21/2 t)作业题解答作业题解答l 第三题第三题 00.511.522.533.5051015