求微分方程的解PPT学习教案

1、会计学1 求微分方程的解求微分方程的解 q 考虑一维经典初值问题 00 , , ( ,)() , dy f x yy xyxa b dx u 基本思想：用差商代替微商 根据 Talyor 公式，y(x) 在点 xk 处有 2 ()()()()() kkk y xy xxxyxOx 1kk hxx 11 ()()()() ( ) kkkk k y xy xy xy xdy O h dxhh x 2 1 ()()()() kkk y xy xhyxO h 第1页/共23页 q 具体步骤： 等距剖分： 0121nn axxxxxb 步长： 1 0 1 21() /, , , kk hxxknnba

2、u 分割求解区间 u 差商代替微商 1 ()() kk k y xy xdy dxh x 1 ()()() kkk y xy xh yx 得方程组：00 1 1 () (,) kkkk kk yy x yyh f xy xxh 分割求解区间，差商代替微商，解代数方程 为分割点 0 n k k x k = 0, 1, 2, ., n-1 yk 是 y (xk) 的近似 第2页/共23页 例：用 Euler 法解初值问题 2 2 0 2 01 , ( ) dyx y dxyx y 取步长 h = (2 - 0)/n = 2/n，得差分方程 00 1 1 0 1 2 , (,)(/) kkkkkkk

3、k kk xy yyh f xyyh yxy xxh 当 h=0.4，即 n=5 时，Matlab 源程序见下一页 解： 第3页/共23页 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 y=y+h*subs(f,x,y,x,y); x=x+h; szj=szj;x,y; end szj plot(szj(:,1),szj(:,2), ) 第4页/共23页 解析解： 1 3 3 52 2 33 / x yex 解析解 近似解 1 3

4、3 52 2 33 / x yex 第5页/共23页 q 为了减小误差，可采用以下方法： u 让步长 h 取得更小一些； u 改用具有较高精度的数值方法： q 龙格-库塔方法 Runge-Kutta (龙格-库塔) 方法 u 是一类求解常微分方程的数值方法 u 有多种不同的迭代格式 第6页/共23页 q 用得较多的是 四阶R-K方法 001 11234 (22)/6 (), kk kk yy xxxh yyh LLLL 1 21 32 43 22 22 (,) (/ ,/ ) (/ ,/ ) (,) kk kk kk kk Lf xy Lf xhyhL Lf xhyhL Lf xh yhL 其

5、中 第7页/共23页 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; end plot(szj(:,1),szj(:,2),

6、 dg-) 第8页/共23页 第9页/共23页 第10页/共23页 q 用 Maltab自带函数 解初值问题 u 求解析解：dsolve （Wang p59） u 求数值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb 第11页/共23页 q dsolve 的使用 y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v) 其中 y 为输出， eq1、eq2、.为微分方程，cond1、cond2、.为初值条件，v 为自变量。 例 1：求微分方程 的通解，并验证。 2 2 x dy xyxe dx y=dsolve(Dy

7、+2*x*y=x*exp(-x2),x) syms x; diff(y)+2*x*y - x*exp(-x2) 第12页/共23页 q 几点说明 l 如果省略初值条件，则表示求通解； l 如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x); dy/dx = 2x dsolve(Dy=2*x); dy/dt = 2x l 若找不到解析解，则返回其积分形式。 l 微分方程中用 D 表示对 自变量 的导数，如： Dy y； D2y y； D3y y 第13页/共23页 例 2：求微分方程 在初值条件 下的特解，并画出解函数的图形。0 x xyye y=dsolve(x*Dy+y-e

8、xp(x)=0,y(1)=2*exp(1),x) ezplot(y); 12( )ye 第14页/共23页 例3：求微分方程组 在初值条件 下的特解，并画出解函数的图形。 5 30 t dx xye dt dy xy dt 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); 0 0 1 0 | | t t x y 注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（structure）类型的数据。 第15页/共23

9、页 例： 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的输出个数只能为一个 或 与方程个数相等。 第16页/共23页 T,Y = solver(odefun,tspan,y0) 其中 y0 为初值条件，tspan为求解区间；Matla

10、b在数值求解时自动对求解区间进行分割，T (向量) 中返回的是分割点的值(自变量)，Y (向量) 中返回的是解函数在这些分割点上的函数值。 solver 为Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、 ode23tb） 没有一种算法可以有效地解决所有的 ODE 问题，因此MATLAB 提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。 第17页/共23页 求解器求解器 ODE类型类型特点特点说明说明 ode45 非刚性非刚性单步法；单步法；4，5 阶阶 R-K 方法；方法； 累计截断误差为累计截断误差为 (x

11、)3 大部分场合的大部分场合的首选方法首选方法 ode23 非刚性非刚性单步法；单步法；2，3 阶阶 R-K 方法；方法； 累计截断误差为累计截断误差为 (x)3 使用于精度较低的情形使用于精度较低的情形 ode113非刚性非刚性多步法；多步法；Adams算法；高低精算法；高低精 度均可到度均可到 10-310-6 计算时间比计算时间比 ode45 短短 ode23t适度刚性适度刚性 采用梯形算法采用梯形算法适度刚性情形适度刚性情形 ode15s刚性刚性多步法；多步法；Gears 反向数值微反向数值微 分；精度中等分；精度中等 若若 ode45 失效时，可失效时，可 尝试使用尝试使用 ode2

12、3s刚性刚性单步法；单步法；2 阶阶Rosebrock 算算 法；低精度法；低精度 当精度较低时，计算时当精度较低时，计算时 间比间比 ode15s 短短 ode23tb刚性刚性梯形算法；低精度梯形算法；低精度当精度较低时，计算时当精度较低时，计算时 间比间比ode15s短短 第18页/共23页 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,

13、0:0.1:0.5,1); %此时 x=0:0.1:0.5 T,Y = solver(odefun,tspan,y0) 求初值问题 的数值解，求解范围为 0,0.5 2 222 01( ) dy yxx dx y 例 4： 第19页/共23页 如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时需用函数文件来定义该常微分方程组。 12 2 2121 12 1 01 00 7 / /() ( ),( ), dxdtx dxdtxxx xx 令 ，则原方程可化为 12 , dy xy x dt 求解 Ver der Pol 初值问题 2 2 2 10 01 00 7 () ( ),( ), d ydy yy dtdt yy 例 5： 第20页/共23页 l 先编写函数文件 verderpol.m function xprime=verderpol(t,x) global mu; xprime=x(2); mu*(1-x(1)2)*x(2) - x(1); l 再编写脚本文件 vdpl.m，在命令窗