第十讲--求微分方程的解

求微分方程的解数学实验自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。本实验主要研究如何用Matlab来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的基本数值解法－－Euler折线法。问题背景和实验目的考虑一维经典初值问题00,,(,)()[,]dyfxyyxyxabdx基本思想：用差商代替微商根据Talyor公式，y(x)在点xk处有2()()()'()()kkkyxyxxxyxOx1kkhxx11()()()()()kkkkkyxyxyxyxdyOhdxhhx21()()'()()kkkyxyxhyxOhEuler折线法初值问题的Euler折线法具体步骤：等距剖分：0121nnaxxxxxb步长：10121()/,,,,,kkhxxknnba分割求解区间差商代替微商1()()kkkyxyxdydxhx1()()'()kkkyxyxhyx得方程组：0011()(,)kkkkkkyyxyyhfxyxxh分割求解区间，差商代替微商，解代数方程为分割点0nkkxk=0,1,2,...,n-1yk是y(xk)的近似Euler折线法举例例：用Euler法解初值问题220201[,]()dyxydxyxy取步长h=(2-0)/n=2/n，得差分方程00211012,(,)(/)kkkkkkkkkkxyyyhfxyyhyxyxxh当h=0.4，即n=5时，Matlab源程序见Euler.m解：Euler折线法源程序clearf=sym('y+2*x/y^2');a=0;b=2;h=0.4;n=(b-a)/h;x=0;y=1;szj=[x,y];fori=1:ny=y+h*subs(f,{'x','y'},{x,y});x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2),'or-')holdonezplot('(5*exp(3*x)/3-2*x-2/3)^(1/3)',[0,2])legend('近似解','解析解')Euler折线法举例（续）解析解：13352233/xyex解析解近似解Runge-Kutta方法为了减小误差，可采用以下方法：让步长h取得更小一些；改用具有较高精度的数值方法：龙格-库塔方法Runge-Kutta(龙格-库塔)方法是一类求解常微分方程的数值方法有多种不同的迭代格式Runge-Kutta方法用得较多的是四阶R-K方法00111234(22)/6(),kkkkyyxxxhyyhLLLL12132432222(,)(/,/)(/,/)(,)kkkkkkkkLfxyLfxhyhLLfxhyhLLfxhyhL其中四阶R-K方法源程序clear;f=sym('y+2*x/y^2');a=0;b=2;h=0.4;n=(b-a)/h;x=0;y=1;szj=[x,y];fori=1:nl1=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];endplot(szj(:,1),szj(:,2),'dg-')Runge-Kutta方法Euler法与R-K法误差比较程序EuRK.mMatlab解初值问题用Maltab自带函数解初值问题求解析解：dsolve求数值解：ode45、ode23、ode113、ode23t、ode15s、ode23s、ode23tbdsolve求解析解dsolve的使用y=dsolve('eq1','eq2',...,'cond1','cond2',...,'v')其中y为输出，eq1、eq2、...为微分方程，cond1、cond2、...为初值条件，v为自变量。例1：求微分方程的通解，并验证。22xdyxyxedxy=dsolve('Dy+2*x*y=x*exp(-x^2)','x')symsx;diff(y)+2*x*y-x*exp(-x^2)dsolve的使用几点说明如果省略初值条件，则表示求通解；如果省略自变量，则默认自变量为tdsolve('Dy=2*x','x');％dy/dx=2xdsolve('Dy=2*x');％dy/dt=2x若找不到解析解，则返回其积分形式。微分方程中用D表示对自变量的导数，如：Dyy'；D2yy''；D3yy'''dsolve举例例2：求微分方程在初值条件下的特解，并画出解函数的图形。0'xxyyey=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y);12()yedsolve举例例3：求微分方程组在初值条件下的特解，并画出解函数的图形。530tdxxyedtdyxydt[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]);0010||ttxy注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（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,options)其中y0为初值条件，tspan为求解区间(自变量的初值和终值)；Matlab在数值求解时自动对求解区间进行分割，T(向量)中返回的是分割点的值(自变量)，Y(向量)中返回的是解函数在这些分割点上的函数值。options设定误差限(缺省:相对10-3,绝对10-6),命令：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分别为设定的相对误差和绝对误差.solver为Matlab的ODE求解器（可以是ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）没有一种算法可以有效地解决所有的ODE问题，因此Matlab提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。ode只能直接求解一阶微分，高阶方法：把一个N阶微分方程转化为N个一阶微分方程组。Matlab提供的ODE求解器求解器ODE类型特点说明ode45非刚性通常情况下最好的单步法；4，5阶R-K方法；误差为(△x)3大部分场合的首选方法(速度慢)ode23非刚性单步法；2，3阶R-K方法；累计截断误差为(△x)3使用于精度较低的情形(比ode45更小步长)ode113非刚性多步法；Adams算法；高低精度均可到10-3～10-6计算时间比ode45短,适合连续的系统ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear’s反向数值微分；精度中等（可用隐式）若ode45失效时，可尝试使用(速度快)ode23s刚性单步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短参数说明odefun为显式常微分方程，可以用命令inline定义，或在函数文件中定义，然后通过函数句柄调用。fun=inline('-2*y+2*x^2+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]222201()dyyxxdxy例4：数值求解举例如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时需用函数文件来定义该常微分方程组。122212112101007//()(),(),dxdtxdxdtxxxxx令，则原方程可化为12,dyxyxdt求解VerderPol初值问题2221001007()(),'(),dydyyydtdtyy例5：数值求解举例先编写函数文件verderpol.mfunctionxprime=verderpol(t,x)%t,x分别给出的时间向量和状态向量globalmu;xprime=[x(2);mu*(1-x(1)^2)*x(2)-x(1)];再编写脚本文件vdpl.m，在命令窗口直接运行该文件。clear;globalmu;mu=7;y0=[1;0];[t,x]=ode45('verderpol',[0,40],y0);plot(t,x(:,1),'rd');holdon[t1,x1]=ode23('verderpol',[0,40],y0);plot(t1,x1(:,1),'bo');先编写函数文件verderpol1.m再编写脚本文件vdpl1.m，在命令窗口直接运行该文件。clear;globalmu;mu=7;y0=[1;0];[t,x]=ode45('verderpol1',[0,40],y0);plot(t,x(:,1),'rd');holdon[t1,x1]=ode15s('verderpol1',[0,40],y0);plot(t1,x1(:,1),'bo');functiondy=verderpol1(t,y)globalmu;dy=zeros(2,1);dy(1)=y(2);dy(2)=mu*(1-y(1)^2)*y(2)-y(1);数值求解举例例6：求解微分方程,1)0(,1'ytyy要求求解析解和数值解，并进行比较。s=dsolve('Dy=-y+t+1','y(0)=1','t')%解析解%数值解，先编写M文件fun1.mfunctionf=fun1(t,y)f=-y+t+1;%m文件中输入clear;t=0:0.1:1;y=t+exp(-t);plot(t,y);holdon%保留已经画好的图形，两个图形合并在一起[t,y]=ode45('fun1',[0,1],1);plot(t,y,'r*')%画数值解图形，用*画xlabel('t')ylabel('y')Matlab求解微分方程小结Matlab函数求解析解（通解或特解），用dsolve求数值解（特解），用ode45、ode23...Matlab编程Euler折线法Runga-Kutta方法例洛伦兹模型由如下常微分方程组描述(蝴蝶效应)