《第八章常微分方程的初值问题课件》由会员分享，可在线阅读，更多相关《第八章常微分方程的初值问题课件（53页珍藏版）》请在金锄头文库上搜索。
1、第八章 常微分方程的初值问题,常微分方程（ODE） ：只含有未知函数的导数的方程。 ODE的阶数： 最高阶导数的阶数,ODE问题,初值问题：,边值问题：,已知初始点处的条件,已知初始点和末端点处的条件,本章只讨论初值问题,一阶ODE问题的形式,1、如果 f 与 y 无关，则计算变为第5章（数值积分）中讨论的任一种直接积分方法应用,初始条件,2、如果 f 是 y 的函数 ，积分过程将不同于前者。,若 f 是 y 的线性函数，如：f=ay+b 其中a,b是常数或是 t 的函数， 此时原方程称为线性ODE,若 f 不是线性函数，方程就称为非线性ODE。,一、求ODE的解析解,dsolve,输出变量列
2、表=dsolve(eq1,eq2, . , eqn, cond1,cond2, . , condn, v1,v2,vn),其中 eq1、eq2、.、eqn为微分方程，cond1、cond2、.、condn为初值条件，v1,v2,vn 为自变量。,注1： 微分方程中用 D 表示对 自变量 的导数，如：,Dy y； D2y y； D3y y,例 ：求微分方程 的通解，并验证。, y=dsolve(Dy+2*x*y=x*exp(-x2),x) 结果为 y =(1/2*x2+C1)*exp(-x2), syms x; diff(y)+2*x*y - x*exp(-x2) 结果为 ans = 0,注2：
3、如果省略初值条件，则表示求通解；,注3：如果省略自变量，则默认自变量为 t,例 ： dsolve(Dy=2*x) %dy/dt = 2x 结果为 ans = 2*x*t+C1,例 ：求微分方程 在初值条件 下的特解，并画出解函数的图形。, y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) 结果为 y = (exp(x)+exp(1)/x ezplot(y) %此时的y已经是符号变量，故不用ezplot(y),dsolve(D3y=-y,y(0)=1,Dy(0)=0,D2y(0)=0) 结果为ans = 1/3*exp(-t)+2/3*exp(1/2*t)*co
4、s(1/2*3(1/2)*t),例：,注4：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（structure）类型的数据。,例：求微分方程组 在初值条件 下的特解,x,y=dsolve(Dx+5*x+y=exp(t), . Dy-x-3*y=0, x(0)=1, y(0)=0, t),x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0, t),或 r=dsolve(Dx+5*x+y=exp(t), Dy-x-3*y=0, . x(0)=1, y(0)=0,
5、t),这里返回的 r 是一个 结构类型 的数据,r.x %查看解函数 x(t) r.y %查看解函数 y(t),dsolve的输出个数只能为一个 或 与方程个数相等。,例 ：用dsolve求解著名的Van der Pol方程, syms mu; y=dsolve(D2y+mu*(y2-1)*Dy+y) 结果：Warning, cannot find an explicit solution y = %N为总的节点数 x(1)=x0;y(1)=y0; for n=1:N-1 x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n); end,下面再用向前欧
6、拉法数值求解，为此先编写向前欧拉法的程序,最后通过图形比较用向前欧拉得到的数值解和解析解的误差,fun=inline(y+2*x/y2,x,y); x1,y1=Euler_bf(fun,0,1,2,0.1); x2,y2=Euler_bf(fun,0,1,2,0.05); x=0:0.01:2; y=1/3*(-18-54*x+45*exp(3*x).(1/3); plot(x1,y1,*,x2,y2,x,x,y) axis(0,2,0,9),向前欧拉方法截断误差为,对,两边从xn 到 xn+1 积分得：,2、改进的Euler法,yn 近似代替y (xn),用梯形代替右边的积分,梯形法,从n=
7、0开始计算，每步都要求解一个关于yn+1的方程（一般是一个非线性方程），可用如下的迭代法计算：,利用此算法，可得：,利用,（ 为允许误差）来控制是否继续迭代,迭代法太麻烦，实际上，当h取得很小时，只让上式中的第二式迭代一次就可以，即,改进的Euler法（也叫欧拉预估校正法）,预估算式,校正算式,改进的Euler法=向前欧拉法+梯形法,向后Euler法依赖于向后差分近似，其格式为：,3、向后Euler法,精度与向前欧拉法相同。如果 f 为非线性函数，则与改进的Euler算法一样，在每一步中需要采用迭代法。该方法有两个优点：（a）绝对稳定；（b）如果解为正，则可保证数值计算结果也为正。,三、龙格库
8、塔（Runge-kutta）方法,Euler法的一个重要缺陷是精度太低。为了获得 高精度，就要减小 h ，这不仅会增加计算时间，也 会产生舍入误差。,1、二阶Runge-kutta方法,或,其实就是欧拉预估校正方法（或改进的欧拉法）,例,用改进的欧拉法（即二阶龙格-库塔法）数值求解并与向前欧拉法、解析解画图比较,前面已求出解析解：,function x,y=Euler_mend(fun,x0,y0,xmax,h) %fun为显式一阶微分方程右端的函数 %x0,y0为初始条件,满足y(x0)=y0 %xmax为x的取值最大值 %h为将x等分后的步长 N=(xmax-x0)/h+1;%N为总的节点
9、数 x(1)=x0;y(1)=y0; for n=1:N-1 x(n+1)=x(n)+h; k1=h*feval(fun,x(n),y(n); k2=h*feval(fun,x(n+1),y(n)+k1); y(n+1)=y(n)+1/2*(k1+k2); end,先编写改进的欧拉法的程序：,再分别调用Euler_bf.m和Euler_mend.m求解：,fun=inline(y+2*x/y2,x,y); x1,y1=Euler_bf(fun,0,1,2,0.1); x2,y2=Euler_mend(fun,0,1,2,0.1); x=0:0.01:2; y=1/3*(-18-54*x+45*
10、exp(3*x).(1/3); plot(x1,y1,*,x2,y2,s,x,y) axis(0,2,0,9),3、三阶龙格库塔方法,常见的三阶龙格库塔方法的格式为：,二阶龙格库塔方法截断误差为,精度不高，希望通过增加计算f的次数提高截断误差的阶数，为此引入,4、四阶龙格库塔方法,常见的四阶龙格-库塔方法有两种， 一种基于1/3辛普森法，格式：,另一种基于3/8辛普森法，格式：,例,用二阶龙格-库塔法和四阶龙格-库塔法数值求解并与、解析解画图比较,前面已求出解析解：,先来编写四阶龙格-库塔法（基于1/3辛普森法）的程序：,function x,y=RK4(fun,x0,y0,xmax,h) %
11、fun为显式一阶微分方程右端的函数 %x0,y0为初始条件,满足y(x0)=y0 %xmax为x的取值最大值 %h为将x等分后的步长 N=(xmax-x0)/h+1;%N为总的节点数 x(1)=x0;y(1)=y0; for n=1:N-1 x(n+1)=x(n)+h; k1=h*feval(fun,x(n),y(n); k2=h*feval(fun,x(n)+1/2*h,y(n)+k1/2); k3=h*feval(fun,x(n)+1/2*h,y(n)+k2/2); k4=h*feval(fun,x(n+1),y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+
12、k4); end,fun=inline(y+2*x/y2,x,y); x1,y1=Euler_mend(fun,0,1,2,0.1); x2,y2=RK4(fun,0,1,2,0.1); x=0:0.1:2; y=1/3*(-18-54*x+45*exp(3*x).(1/3); plot(x1,y1,*,x2,y2,s,x,y) axis(0,2,0,9),再分别调用Euler_mend.m和RK4.m求解：,从图形上看，似乎二阶龙格库塔法与四阶龙格-库塔法同样接近解析解，故再从数值结果比较看看哪种误差小,err=abs(y-y1),abs(y-y2),结果为 err = 0 0 0.0009
13、 0.0000 0.0016 0.0000 0.0022 0.0000 0.0026 0.0000 0.0031 0.0000 0.0035 0.0000 0.0040 0.0000 0.0047 0.0000 0.0054 0.0000 0.0062 0.0000 0.0073 0.0000 0.0084 0.0000 0.0098 0.0000 0.0115 0.0000 0.0133 0.0000 0.0155 0.0000 0.0181 0.0000 0.0210 0.0000 0.0243 0.0000 0.0281 0.0000,四、二阶ODE问题,二阶ODE的一般形式为：,其中
14、是常数或是 的函数，后两个方程为初 始条件。,如果 与u无关，则该方程称为线性ODE;,如果三者之中有一个是u或 的函数，称为非线性的。,下面着重研究用向前Euler法求解二阶ODE，及MATLAB程序。,所以二阶ODE分解为两个一阶ODE：,设：,定义,则上述两个一阶ODE写为：,其向前Euler法的格式为：,例 求解二阶ODE,解：设,令,则原方程的向量形式为：,向前Euler法的格式为：,h = 0.05; t_max=5; n=1; y(:,1) = 0; 1;t(1) = 0; while t(n) t_max y(:,n+1) = y(:,n) + h*f_def(y(:,n),t
15、); t(n+1) = t(n)+h ; n=n+1; end axis(0 5 -1 1);plot(t, y(1,:), t, y(2,:),:),function f=f_def(y,t) f=y(2);(-5*abs(y(2)*y(2)-20*y(1);,先用函数文件定义f(u,v,t),再用向前Euler法求解,五、matlab相关命令数值求解ODE,X,Y = 求解函数(fun,x0,xf,y0,option，p1,p2,.),其中：（1）fun是用inline或函数文件定义的显式常微分方程的函数名,函数文件格式：,或inline格式：,function yd=函数名(x,y,flag,p1,p2,) yd=.,函数名=inline(显式方程,x,y,flag,p1,p2,.),x为自变量，y为因变量，yd为y的导数, flag用于控制求解过程，p1，p2为方程中的参数,X,Y = 求解函数(fun,x0,xf,y0,option，p1,p2,.),其中： （2）x0,xf为求解区间， y0 为初值条件； （3）option（可省略）为控制选项（用于设置误差限，求解方程最大允许步长等等）， （4）p1,p2为微分方程中的附加参数 （5）Matlab在数值求解时自动对求
