《《求微分方程的解》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《求微分方程的解》PPT课件.ppt（34页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、实验四求微分方程的解数学实验数学实验1q 自牛顿发明微积分以来，微分方程在描述事物运自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。数学建模所得到的方程，绝大多数是微分方程。q 由于实际应用的需要，人们必须求解微分方程。由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。数微分方程需要利用数值方法来近似求解。q 本实验主要研究如何用本实验主要研究如何用 Ma
2、tlab 来计算微分方程来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的（组）的数值解，并重点介绍一个求解微分方程的基本数值解法基本数值解法Euler折线法折线法。问题背景和实验目的问题背景和实验目的2q 考虑一维经典考虑一维经典初值问题初值问题u 基本思想：基本思想：用差商代替微商用差商代替微商根据根据 Talyor 公式，公式，y(x)在点在点 xk 处有处有Euler 折线法折线法3初值问题的初值问题的Euler折线法折线法q 具体步骤：具体步骤：等距剖分：等距剖分：步长：步长：u 分割求解区间分割求解区间分割求解区间，差商代替微商，解代数方程分割求解区间，差商代替微商，解代数
3、方程 为分割点为分割点u 差商代替微商差商代替微商得迭代格式：得迭代格式：k=0,1,2,.,n-1yk 是是 y(xk)的近似的近似4Euler 折线法举例折线法举例例：例：用用 Euler 法解初值问题法解初值问题取步长取步长 h=(2-0)/n=2/n，得差分方程得差分方程当当 h=0.4，即即 n=5 时，时，Matlab 源程序见源程序见 fuluA.m解：解：5Euler 折线法源程序折线法源程序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=
4、y+h*subs(f,x,y,x,y);x=x+h;szj=szj;x,y;endszjplot(szj(:,1),szj(:,2),or-)6 Euler折线法举例（续）折线法举例（续）解析解：解析解：解析解解析解近似解近似解y=1/3*(-18-54*x+45*exp(3*x)(1/3)7Runge-Kutta 方法方法q 为了减小误差，可采用以下方法：为了减小误差，可采用以下方法：u 让步长让步长 h 取得更小一些；取得更小一些；u 改用具有较高精度的数值方法：改用具有较高精度的数值方法：q 龙格龙格-库塔方法库塔方法Runge-Kutta(龙格龙格-库塔库塔)方法方法u 是是一类一类求
5、解常微分方程的数值方法求解常微分方程的数值方法u 有多种不同的迭代格式有多种不同的迭代格式8Runge-Kutta 方法方法q 用得较多的是用得较多的是 四阶四阶R-K方法方法其中其中9四阶四阶 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+
6、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-)10Runge-Kutta 方法方法11Euler 法与法与 R-K法误差比较法误差比较12Matlab 解初值问题解初值问题q 用用 Maltab自带函数自带函数 解初值问题解初值问题u 求解析解：求解析解：dsolveu 求数值解：求数值解：ode45、ode23、ode113、ode23t、ode15s、ode23s、ode23tb13dsolve 求解析解求解析解q dsolve 的使用的使用y=dsolve(eq1,e
7、q2,.,cond1,cond2,.,v)其中其中 y 为输出，为输出，eq1、eq2、.为微分方程，为微分方程，cond1、cond2、.为初值条件，为初值条件，v 为自变量。为自变量。例例 1：求微分方程求微分方程 的通解，并验证。的通解，并验证。y=dsolve(Dy+2*x*y=x*exp(-x2),x)syms x;diff(y)+2*x*y-x*exp(-x2)14dsolve 的使用的使用q 几点说明几点说明l 如果省略初值条件，则表示求通解；如果省略初值条件，则表示求通解；l 如果省略自变量，则默认自变量为如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x);d
8、y/dx=2xdsolve(Dy=2*x);dy/dt=2xl 若找不到解析解，则返回其积分形式。若找不到解析解，则返回其积分形式。l 微分方程中用微分方程中用 D 表示对表示对 自变量自变量 的导数，如：的导数，如：Dy y；D2y y；D3y y15dsolve 举例举例例例 2：求微分方程求微分方程 在初值条件在初值条件 下的特解，并画出解函数的图形。下的特解，并画出解函数的图形。y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x)ezplot(y);16dsolve 举例举例例例3：求微分方程组求微分方程组 在初值条件在初值条件 下的特解，并画出解函数的图
9、形。下的特解，并画出解函数的图形。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);注：解微分方程组时，如果所给的输出个数与方程个数相同，注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解则方程组的解按词典顺序按词典顺序输出；如果只给一个输出，则输出输出；如果只给一个输出，则输出的是一个包含解的的是一个包含解的结构结构（structure）类型的数据。类型的数据。17 dsolve 举例举例例：例：x,y=dsolve(Dx+5*x=0,Dy-3*y=0,.x(0)=1,y(0)=1
10、,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的输出个数只能为一个的输出个数只能为一个 或或 与方程个数相等与方程个数相等18数值求解数值求解T,Y=solver(odefun,tspan,y0)其中其中 y0 为初值
11、条件，为初值条件，tspan为求解区间；为求解区间；Matlab在数值求解在数值求解时时自动对求解区间进行分割自动对求解区间进行分割，T(列向量列向量)中返回的是分割点中返回的是分割点的值的值(自变量自变量)，Y(数组数组)中返回的是这些分割点上的近似解，中返回的是这些分割点上的近似解，其列数等于因变量的个数。其列数等于因变量的个数。solver 为为Matlab的的ODE求解器求解器（可以是（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）没有一种算法可以有效地解决所有的没有一种算法可以有效地解决所有的 ODE 问题，因此问题，因此MA
12、TLAB 提供了多种提供了多种ODE求解器求解器，对于不同的对于不同的ODE，可以调用不同的可以调用不同的求解器求解器。19Matlab的的ODE求解器求解器求解器求解器 ODE类型类型特点特点说明说明ode45非刚性非刚性单步法；单步法；4 4，5 5 阶阶 R-K R-K 方法；方法；累计截断误差为累计截断误差为 (x)3大部分场合的大部分场合的首选方法首选方法ode23非刚性非刚性单步法；单步法；2 2，3 3 阶阶 R-K R-K 方法；方法；累计截断误差为累计截断误差为 (x)3使用于精度较低的情形使用于精度较低的情形ode113非刚性非刚性多步法；多步法；Adams算法；高低精算法
13、；高低精度均可到度均可到 10-310-6计算时间比计算时间比 ode45 短短ode23t适度刚性适度刚性 采用梯形算法采用梯形算法适度刚性情形适度刚性情形ode15s刚性刚性多步法；多步法；Gears Gears 反向数值微反向数值微分；精度中等分；精度中等若若 ode45 失效时，可失效时，可尝试使用尝试使用ode23s刚性刚性单步法；单步法；2 2 阶阶Rosebrock Rosebrock 算法；算法；低精度低精度当精度较低时，计算时当精度较低时，计算时间比间比 ode15s 短短ode23tb刚性刚性梯形算法；低精度梯形算法；低精度当精度较低时，计算时当精度较低时，计算时间比间比o
14、de15s短短20参数说明参数说明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=solver(odefun,tspan,y0)求求 的数值解，求解区间为的数值解，求解区
15、间为 0,0.5例例：21数值求解举例数值求解举例如果需求解的问题是如果需求解的问题是高阶高阶常微分方程，则需将其化为常微分方程，则需将其化为一阶常一阶常微分方程组微分方程组，此时必须用，此时必须用函数文件函数文件来定义该常微分方程组。来定义该常微分方程组。令令 求解求解 Ver der Pol 初值问题初值问题例例：22数值求解举例数值求解举例l 先编写函数文件先编写函数文件 verderpol.mfunction xprime=verderpol(t,x)global mu;xprime=x(2);mu*(1-x(1)2)*x(2)-x(1);l 再编写脚本文件再编写脚本文件 vdpl.m
16、，在命令窗口直接运行该文件，在命令窗口直接运行该文件 clear;global mu;mu=7;y0=1;0;%t,x=ode45(verderpol,0,40,y0);t,x=ode45(verderpol,0,40,y0);plot(t,x(:,1);23求解微分方程小结求解微分方程小结q Matlab 函数函数u 求解析解（求解析解（通解通解或特解），用或特解），用 dsolveu 求数值解（特解），用求数值解（特解），用 ode45、ode23.q Matlab 编程编程u Euler 折线法折线法u Runga-Kutta 方法方法24如何定义函数文件如何定义函数文件T,Y=ode4
17、5(odefun,tspan,y0)T,Y=ode23(odefun,tspan,y0)当当 是函数向是函数向量时呢？量时呢？微分方程组微分方程组odefun那么那么odefun就是就是25如何定义函数文件如何定义函数文件T,Y=ode45(odefun,tspan,y0)T,Y=ode23(odefun,tspan,y0)ode45、ode23 等函数可用于求解等函数可用于求解显式常微分方程显式常微分方程当当 是向量函数时，所对应的方程即为是向量函数时，所对应的方程即为微分方程组微分方程组odefun26举例说明举例说明fun=inline(-2*y+2*x2+2*x,x,y);x,y=od
18、e23(fun,0,0.5,1);例：求初值问题例：求初值问题 的数值解。的数值解。解法一：解法一：使用使用 inline 定义微分方程定义微分方程 odefunodefun 为方程右端项为方程右端项 f(t,y)可以用可以用 inline 定义（只适合于定义（只适合于单个方程单个方程的情形的情形）通过通过函数文件函数文件定义，然后用函数句柄调用（适合所有情形）定义，然后用函数句柄调用（适合所有情形）注：自变量必须在前面，因变量在后面！注：自变量必须在前面，因变量在后面！27举例说明（单个方程）举例说明（单个方程）function dy=myfun1(x,y)dy=-2*y+2*x2+2*x;
19、解法二：解法二：通过通过函数文件函数文件定义微分方程定义微分方程 odefun1、先编写函数文件、先编写函数文件 myfun1.mclear;x,y=ode23(myfun1,0,0.5,1);2、编写主文件编写主文件 main1.m或直接在或直接在 Matlab 命令窗口输入上面的语句。命令窗口输入上面的语句。28举例说明（方程组）举例说明（方程组）解解：此时只能通过：此时只能通过函数文件函数文件定义微分方程定义微分方程 odefun例：求例：求 ,的数值解。的数值解。function dy=myfun2(t,y)dy=zeros(3,1);%dy must be a column vect
20、or!dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);1、先编写函数文件、先编写函数文件 myfun2.mclear;T,Y=ode45(myfun2,0,12,0,1,1);2、编写主文件编写主文件 main2.mdy=y(2)*y(3);-y(1)*y(3);.-0.51*y(1)*y(2);29思考思考function dy=myfun2(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);1、函数文件、函数文件 myfun2.m能不
21、能写成下面形式？能不能写成下面形式？function dy=myfun2(t,x,y,z)dy=zeros(3,1);dy(1)=y*z;dy(2)=-x*z;dy(3)=-0.51*x*y;X30说明说明odefun变量属性必须变量属性必须一一对应！一一对应！function dy=myfun2(t,y)如果是常微分方程组，如果是常微分方程组，y 就是列向量！就是列向量！dy 必须是必须是列列向量向量，长度为方程组的个数，通常与，长度为方程组的个数，通常与y的长度相同！的长度相同！函数中的输入参数和输出参数是函数中的输入参数和输出参数是形参形参，名字可以任意取，但必须，名字可以任意取，但必须
22、满足上述条件。即输入参数有两个，第一个表示满足上述条件。即输入参数有两个，第一个表示自变量自变量，第二个，第二个是由是由因变量因变量组成的组成的列向量列向量，输出参数必须是列向量。，输出参数必须是列向量。31例例function dy=myfun3(t,y)dy=zeros(2,1);dy(1)=y(2)+t;dy(2)=t-2;例：解初值问题：例：解初值问题：，function out=myfun3(t,y)out=y(2)+t;t 2;function yprime=myfun3(x,y)yprime=y(2)+x;x 2;clear;T,Y=ode45(myfun3,0,10,1,1);
23、2、主文件主文件 main3.m1、函数文件、函数文件 myfun3.m32高阶常微分方程高阶常微分方程高阶高阶常微分方程常微分方程例：例：Van der Pol 初值问题初值问题令令 ，则原方程可化为，则原方程可化为一阶一阶常微分方程组常微分方程组变量替换变量替换化为化为参数怎么处理？参数怎么处理？用全局变量传递用全局变量传递33参数传递参数传递1、函数文件、函数文件 verderpol.mfunction xprime=verderpol(t,x)global mu;xprime=x(2);mu*(1-x(1)2)*x(2)-x(1);2、主文件、主文件 vdp1.mclear;global mu;y0=1;0;mu=7;t,x=ode45(verderpol,0,40,y0);%t,x=ode45(verderpol,0,40,y0);plot(t,x(:,1),r-,t,x(:,2),b-);34