计算机应用基础-5-常微分方程求解.ppt

1、第五章 常微分方程的求解,5.1 常微分方程基本特征 5.2 常微分方程的初值问题求解 5.3 常微分方程的边值问题求解 5.4 高阶常微分方程的求解,5.1 常微分方程基本特征,常微分方程（组）包括初值问题和边值问题两种： （1） 初值问题： 初始条件已知，形式如下,对于初值问题方程组，由多个常微分方程及其初始条件构成,（2） 边值问题：已知首末两个端点值，形式如下：,边值问题求解： 打靶法：将边值问题转化成初值问题 松弛法：包括有限差分法和配置法,5.1 常微分方程基本特征,一 欧拉公式,对于常微分方程初值问题,在求解区间a,b上作等距分割，步长 记xn = xn-1+h , n =1,2

2、,m。用差商近似导数计算常微分方程。 做 的在x = x0处的一阶向前差商得： 于是得到：,5.2 常微分方程的初值问题求解,故y(x1)的近似值y1可按 求得。类似地，由 得到计算近似值的向前欧拉公式,5.2 常微分方程的初值问题求解,由yn直接算出yn+1值的计算格式称为显式格式，向前欧拉公式是显式格式。而欧拉方法的几何意义是以y1作为斜率，通过点(x0, y0)做一条直线，它与直线x=x1的交点就是y1。依此类推，是以yn+1作为斜率，经过点的直线与直线x=xn+1的交点。故欧拉法也称为欧拉折线法，如下图所示。,5.2 常微分方程的初值问题求解,龙格-库塔法是求解常微分方程较常用的一种方

3、法，它通过巧妙的线性组合，在显式格式的情况下获得理想的计算精度，大大提高了计算速度。该方法的推导过程不是本课程要研究的对象，本课程主要研究其实际应用,，故直接给出各类龙格-库塔公式。 1、 二阶龙格-库塔,二 龙格-库塔方法,5.2 常微分方程的初值问题求解,二阶龙格-库塔公式的局部截断误差为O(h3), 是二阶精度的计算公式。也可建立高阶的龙格-库塔公式，如三阶龙格-库塔公式、四阶龙格-库塔公式。较常用的是四阶龙格-库塔公式，四阶龙格-库塔公式的局部截断误差界为O(h5), 是四阶精度的计算公式。,5.2 常微分方程的初值问题求解,5.2 常微分方程的初值问题求解,2 三阶龙格-库塔,3、四

4、阶龙格库塔公式,5.2 常微分方程的初值问题求解,龙格-库塔方法的调用格式：,t, y=ode45(ODEfun, tspan, y0, options, p1),自定义的函数,求解器,积分限或离散的向量,初值,Options：积分参数值，不选为默认值 p1: 传递给ODEfun的参数,输出： t: 时间变量，因变量，列向量 y: 状态变量，求解,5.2 常微分方程的初值问题求解,例4-1 求解初值问题：,计算结果： ode45 ode23 精确值 x=1 1.73205081676653 1.73215487941780 1.7320508,5.2 常微分方程的初值问题求解,求解程序”xOD

5、E.m”,龙格-库塔方法步长的选择,怎样选取合适的步长，这在实际计算中是很重要的。由于步长愈小，每步计算的截断误差就愈小；但在一定的求解范围内，需要完成的步数就愈多，这不但引起计算量的增加，而且计算步数的增加往往造成舍入误差的严重积累，反而降低了计算精度。上面介绍的龙格-库塔方法是对定步长（即步长h为常数）而言的，但为了保证精确度，一种有效的措施是在计算过程中自动进行步长调整，即变步长技巧。,5.2 常微分方程的初值问题求解,下面以四阶龙格-库塔方法为例，说明如何自动选择步长，使计算结果满足给定精度的要求。 设从节点xn出发，先以h为步长，利用四阶龙格-库塔公式方法经过一步计算得y(xn+1)

6、的近似值，记为 ，由于公式的局部截断误差是y(h5)，故有 当h不大时，c可近似地看作常数。然后将步长h对折，即取h/2为步长，从出发经过两步计算求y(xn+1)的近似值，记为 ，每一步计算的局部截断误差为c(h/2)5，于是就有,5.2 常微分方程的初值问题求解,把它与上式相比，可得 经整理可得 这表明以 作为y(xn+1) 的近似值，其误差可用先后两次计算结果之差来表示，因而，只需考察 是否成立。若成立，则可将 作为y(xn+1)的近似值；若不成立，则将步长再次对折进行计算，直到不等式成立为止，并取最后的 作为计算结果。以上方法就是计算过程中自动选择步长的方法，也称为变步长方法。从表面上看

7、，为了选择适当的步长，每一步的计算量增加了，但从总体考虑，工作量往往还是减少的。,5.2 常微分方程的初值问题求解,龙格-库塔方法的主要优点是计算精确度较高，能满足通常的计算要求；容易编制程序；每次计算y(xn+1) ，只用到前一步的计算结果yn，因此，在已知初始值y0的条件下，就可自动地进行计算，是单步法，而且计算过程中可随时改变步长。它的缺点是每前进一步需要计算多次f(x,y)的值，因此，计算工作量较大，且其截断误差难以估计。 在实际应用上，一般当要求更高精确度时，采用的办法是缩小步长，而不是采用更高阶的公式，因为高阶公式的计算太复杂，一般选用标准四阶龙格-库塔方法即可。,5.2 常微分方

8、程的初值问题求解,下面用一个例子说明： 用标准四阶龙格-库塔方法计算得 由积分法算得y(2)=2.23607,当h=0.5时相差0.00013，而h=0.25误差小于0.000001，但当 h=0.125时虽然足够精确但计算次数却比h=0.25多了一倍。因此合理选择步长既能保证精度又能减少计算量。,5.2 常微分方程的初值问题求解,5.2 常微分方程的初值问题求解,在区间【0,12】求解下列微分方程组的解,ME_5_13.m,【例5-1】,5.2 常微分方程的初值问题求解,ME_5_5.m,求解Lorenz模型的状态方程 初值为x1(0)=x2(0)=0, x3(0)=10-10,边值问题求解

9、之一-打靶算法,5.3 常微分方程的边值问题求解,二阶微分方程的边值问题的 求解： 在区间a;b上求解该方程的 解，且在这两个边界点上满足： 上式即为边界条件,5.3 常微分方程的边值问题求解,对于线性微分方程 其中 p(x), q(x),f(x)均为给定的函数，边界条件 可简化为： y(a)=a y(b)=b,5.3 常微分方程的边值问题求解,打靶法的步骤： （1）求出下面方程初值问题的数值解 y1(b) (2)求出下面方程初值问题的数值解 y1(b) (3)求出下面方程初值问题的数值解 yp(b),（4） 若y2(b) =0,则计算,5.3 常微分方程的边值问题求解,(5) 计算下面初值问

10、题的 数值解，则y(x)即为原边值问题的数值解,得出一阶微分方程组模型： 取变量 x1=y, x2=y 得到：,function t, y=shooting(f1, f2, tspan, x0f, varargin) t0=tspan(1); tfinal=tspan(2); ga=x0f(1); gb=x0f(2); t, y1=ode45(f1, tspan, 1; 0, varargin); t, y2=ode45(f1, tspan, 0; 1, varargin); t, yp=ode45(f2, tspan, 0; 0, varargin); m=(gb-ga*y1(end, 1)

11、-yp(end,1)/y2(end, 1); t, y=ode45(f2, tspan, ga; m, varargin);,将上述文件存为“shooting.m”,5.3 常微分方程的边值问题求解,编写求解通用函数,5.3 常微分方程的边值问题求解,其中： tspan为初始和终止计算时间构成的向量 为边界值,辅助函数f1( ),辅助函数f2( ),【例5-2】,function xdot=c7fun1(t, x) xdot=x(2); -2*x(1)+3*x(2);,function xdot=c7fun2(t, x) xdot=x(2); t-2*x(1)+3*x(2);,编写两个函数，分

12、别存为m文件,5.3 常微分方程的边值问题求解,C7fun1.m,C7fun2.m,求解：, t, y=shooting(c7fun1, c7fun2, 0, 1, 1, 2) ; plot (t,y) ME_5_10,5.3 常微分方程的边值问题求解,边值问题求解之二-有限差分算法,5.3 常微分方程的边值问题求解,有限差分法是用差分格式近似替代ODE方程中的导数项，从而将ODE离散为线性或非线性方程组而求解，差分格式在偏微分方程中讲解。,5.3 常微分方程的边值问题求解,其中n为中间计算点的个数 i, i=1,2,n-1 为微分方程的数值解 yi+1= i, y1=a和yn=b 其它参数的

13、计算,5.3 常微分方程的边值问题求解,Function x,y=fdiff(funcs, tspan, x0f,n) t0=tspan(1); tfinal=tspan(2); ga=x0f(1); gb=x0f(2); for i=1:n x(i)=t0+h*(i-1); end x0=x(1:n-1); t=-2+h2*feval(funcs,x0,2); tmp= feval(funcs,x0,1); v=1+h*tem/2; w=1-h*tem/2; b= h2*feval(funcs,x0,3);,b(1)=b(1)-w(1)*ga; b(n-1)=b(n-1)-v(n-1)*gb

14、; b=b; A=diag(t); for i=1:n-2 A(i,i+1)=v(i); A(i+1,i)=w(i+1); end y=inv(A)*b; x=x, tfinal; y=ga; y; gb;,5.4 微分方程的解析解方法,边值问题求解之三-符号算法,y=dsolve(f1,f2,fm) y=dsolve(f1,f2,fm,x) % 指明自变量,求 x(t)=x(t)(1-x2(t)的解析解,Syms x t X=dsolve(Dx=x*(1-x2),例5-8,X = 1/(1+exp(-2*t)*C1)(1/2) -1/(1+exp(-2*t)*C1)(1/2),【例5-5】输

15、入信号,syms t y u=exp(-5*t)*cos(2*t+1)+5; uu=5*diff(u,t,2)+4*diff(u,t)+2*u,uu =87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10,5.4 微分方程的解析解方法,求以下微分方程的通解：,步骤1,步骤2,Syms t y y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=, . 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10),y = -87/52*cos(1)*cos(t)2*exp(-5*t)

16、-87/260*cos(1)*sin(t)*cos(t)*exp(-5*t)+87/104*exp(-5*t)*cos(1)+23/26*exp(-5*t)*sin(1)-87/520*sin(1)*cos(2*t)*exp(-5*t)+23/130*cos(1)*cos(2*t)*exp(-5*t)-23/65*sin(1)*sin(t)*cos(t)*exp(-5*t)+5/12+87/104*sin(1)*sin(2*t)*exp(-5*t)-23/26*cos(1)*sin(2*t)*exp(-5*t)-23/13*sin(1)*cos(t)2*exp(-5*t)+C1*exp(-4*

17、t)+C2*exp(-3*t)+C3*exp(-2*t)+C4*exp(-t),5.4 微分方程的解析解方法,syms t y y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=, . 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10, y(0)=3, . Dy(0)=2, D2y(0)=0, D3y(0)=0),y(t)= -87/52*cos(1)*cos(t)2*exp(-5*t)-87/520*sin(1)*cos(2*t)*exp(-5*t)+23/130*cos(1)*cos(2*t)*exp(-5*t)

18、-23/65*sin(1)*sin(t)*cos(t)*exp(-5*t)-87/260*cos(1)*sin(t)*cos(t)*exp(-5*t)+87/104*sin(1)*sin(2*t)*exp(-5*t)-23/26*cos(1)*sin(2*t)*exp(-5*t)-23/13*sin(1)*cos(t)2*exp(-5*t)+5/12+87/104*exp(-5*t)*cos(1)+23/26*exp(-5*t)*sin(1)+(-271/30*cos(1)+41/15*sin(1)-25/4)*exp(-4*t)+(179/8*cos(1)+5/8*sin(1)+73/3)*

19、exp(-3*t)+(-445/26*cos(1)-51/13*sin(1)-69/2)*exp(-2*t)+(133/30*cos(1)+97/60*sin(1)+19)*exp(-t),5.4 微分方程的解析解方法,假设 y(0)=3, y(0)=2, y(0)=y(0)=0,x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t), Dy=4*x+3*y+4*exp(-t),x = -1/6*(-24*C1*exp(-t)+6*C1*exp(1+6(1/2)*t)-3*C1*6(1/2)*exp(1+6(1/2)*t)+6*C1*exp(-(-1+6(1/2)*t)+3*C1*

20、6(1/2)*exp(-(-1+6(1/2)*t)+5*C2*6(1/2)*exp(-(-1+6(1/2)*t)+12*C2*exp(-(-1+6(1/2)*t)-5*C2*6(1/2)*exp(1+6(1/2)*t)+12*C2*exp(1+6(1/2)*t)-24*C2*exp(-t)-2*C3*6(1/2)*exp(-(-1+6(1/2)*t)-6*C3*exp(-(-1+6(1/2)*t)+2*C3*6(1/2)*exp(1+6(1/2)*t)-6*C3*exp(1+6(1/2)*t)+12*C3*exp(-t)-150*exp(-t)+72*exp(-t)*t)/(-2+6(1/2)

21、/(2+6(1/2),5.4 微分方程的解析解方法,y = 2/3*(3*C1*exp(1+6(1/2)*t)-6*C1*exp(-t)+3*C1*exp(-(-1+6(1/2)*t)+C2*6(1/2)*exp(-(-1+6(1/2)*t)+3*C2*exp(-(-1+6(1/2)*t)-C2*6(1/2)*exp(1+6(1/2)*t)+3*C2*exp(1+6(1/2)*t)-6*C2*exp(-t)+3*C3*exp(-t)-C3*6(1/2)*exp(-(-1+6(1/2)*t)+C3*6(1/2)*exp(1+6(1/2)*t)+18*exp(-t)*t-36*exp(-t)/(-

22、2+6(1/2)/(2+6(1/2),5.4 微分方程的解析解方法,5.5 高阶常微分方程求解,一 单个高阶常微分方程处理方法,一个高阶微分方程的一般形式为： 输出变量y(t)的各阶导数初始值为： 选出一组状态变量：,5.5高阶常微分方程求解,原高阶常微分方程变换为一阶微分方程组 初值为：,【例5-3】,5.5 高阶常微分方程求解,function y=vdp_eq(t, x, flag, mu) y=x(2) ; -mu*(x(1)2-1)*x(2)-x(1);,建立函数“vdp_eq.m”,x0=-0.2; -0.7; t_final=20; mu=1; t1, y1=ode45(vdp_eq, 0, t_final, x0, , mu); mu=2; t2, y2=ode45(vdp_eq, 0, t_final, x0, , mu); plot(t1, y1, t2, y2, :) figure; plot(y1(:, 1), y1(:, 2), y2(:, 1), y2(:, 2), :),已知：y(0)=-0.2, y(0)=-0.7, 用数值方法求 Van der pol 方程的解：,ME_5_16.m,5.5 高阶常微分方程求解,5.5 高阶常微分方程求解,二 高阶常微分方程组,5.5 高阶