资源描述
,单击此处编辑母版标题样式,单击此处编辑母版文本样式,*,求微分方程的解,自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。,由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。,本实验主要探讨如何用,Matlab,来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的基本数值解法,Euler,折线法,。,问题背景和实验目的,考虑一维经典,初值问题,基本思想：,用差商代替微商,根据,Talyor,公式，,y,(,x,),在点,x,k,处有,Euler,折线法,初值问题的,Euler,折线法,具体步骤：,等距剖分：,步长：,分割求解区间,差商代替微商,得方程组：,分割求解区间，差商代替微商，解代数方程,为分割点,k,=0,1,2,.,n,-1,y,k,是,y,(,x,k,),的近似,Euler,折线法举例,例：,用,Euler,法解初值问题,取步长,h,=(2-0),/,n,=2,/,n,，,得差分方程,当,h=0.4,，,即,n=5,时，,Matlab,源程序见,fulu1.m,解：,Euler,折线法源程序,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),or-,),Euler,折线法举例（续）,解析解：,解析解,近似解,Runge-Kutta,方法,为了减小误差，可采用以下方法：,让步长,h,取得更小一些；,改用具有较高精度的数值方法：,龙格,-,库塔方法,Runge-Kutta(,龙格,-,库塔,),方法,是,一类,求解常微分方程的数值方法,有多种不同的迭代格式,Runge-Kutta,方法,用得较多的是,四阶,R-K,方法,（教材第,92,页）,其中,四阶,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;,end,plot(szj(:,1),szj(:,2),dg-,),Runge-Kutta,方法,Euler,法与,R-K,法误差比较,Matlab,解初值问题,用,Maltab,自带函数,解初值问题,求解析解：,dsolve,求数值解：,ode45、ode23、ode113、ode23t,、,ode15s、,ode23s、ode23tb,dsolve,求解析解,dsolve,的使用,y=dsolve(,eq1,eq2,.,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),dsolve,的使用,几点说明,如果省略初值条件，则表示求通解；,如果省略自变量，则默认自变量为,t,dsolve(,Dy=2*x,x,);,dy/dx=2x,dsolve(,Dy=2*x,);,dy/dt=2x,若找不到解析解，则返回其积分形式。,微分方程中用,D,表示对,自变量,的导数，如：,Dy y,；,D2y y,；,D3y y,dsolve,举例,例 2：,求微分方程 在初值条件 下的特解，并画出解函数的图形。,y=dsolve(,x*Dy+y-exp(x)=0,y(1)=2*exp(1),x,),ezplot(y);,dsolve,举例,例3：,求微分方程组,在初值条件,下的特解，并画出解函数的图形。,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）,类型的数据。,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,的输出个数只能为一个 或 与方程个数相等。,Matlab,函数数值求解,T,Y=,solver,(,odefun,tspan,y0,),其中,y0,为初值条件，,tspan,为求解区间；,Matlab,在数值求解时,自动对求解区间进行分割,，,T,(,向量,),中返回的是分割点的值,(,自变量,),，,Y,(,向量,),中返回的是解函数在这些分割点上的函数值。,solver,为,Matlab,的,ODE,求解器,（可以是,ode45,、,ode23,、,ode113,、,ode15s,、,ode23s,、,ode23t,、,ode23tb,）,没有一种算法可以有效地解决所有的,ODE,问题，因此,MATLAB,提供了多种,ODE,求解器,，,对于不同的,ODE,，可以调用不同的,求解器,。,Matlab,提供的,ODE,求解器,求解器,ODE,类型,特点,说明,ode45,非刚性,单步法；4，5 阶,R-K,方法；累计截断误差为(,x),3,大部分场合的,首选方法,ode23,非刚性,单步法；2，3 阶,R-K,方法；累计截断误差为(,x),3,使用于精度较低的情形,ode113,非刚性,多步法；,Adams,算法；高低精度均可到 10,-3,10,-6,计算时间比,ode45,短,ode23t,适度刚性,采用梯形算法,适度刚性情形,ode15s,刚性,多步法；,Gears,反向数值微分；精度中等,若,ode45,失效时，可尝试使用,ode23s,刚性,单步法；2 阶,Rosebrock,算法；低精度,当精度较低时，计算时间比,ode15s,短,ode23tb,刚性,梯形算法；低精度,当精度较低时，计算时间比,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.5,T,Y=,solver,(,odefun,tspan,y0,),求初值问题,的数值解，求解范围为 0,0.5,例 4：,数值求解举例,如果需求解的问题是,高阶,常微分方程，则需将其化为,一阶常微分方程组,，此时需用,函数文件,来定义该常微分方程组。,令,，则原方程可化为,求解,Ver der Pol,初值问题,例 5：,数值求解举例,先编写函数文件,verderpol.m,function xprime=verderpol(t,x),global mu;,xprime=x(2);mu*(1-x(1)2)*x(2)-x(1);,再编写脚本文件,vdpl.m,，在命令窗口直接运行该文件。,clear;,global mu;mu=7;,y0=1;0;,t,x=ode45(,verderpol,0,40,y0,),;,plot(,t,x(:,1),r-,),;,Matlab,求解微分方程小结,Matlab,函数,求解析解（,通解,或特解），用,dsolve,求数值解（特解），用,ode45,、,ode23,.,Matlab,编程,Euler,折线法,Runga-Kutta,方法,教材,P97：,练习 1、2、3、4、5、6、7,上机作业,要求：,请在上机之前将程序写好（参考附录），上机时直接输入或修改附录程序即可。,
展开阅读全文