下载

1下载券

加入VIP
  • 专属下载特权
  • 现金文档折扣购买
  • VIP免费专区
  • 千万文档免费下载

上传资料

关闭

关闭

关闭

封号提示

内容

首页 求微分方程的解-2013050914315473

求微分方程的解-2013050914315473.ppt

求微分方程的解-2013050914315473

152*****014@sina.cn
2013-07-23 0人阅读 举报 0 0 0 暂无简介

简介:本文档为《求微分方程的解-2013050914315473ppt》,可适用于人文社科领域

一、问题背景和实验目的二、求解微分方程的若干计算方法三、求解微分方程的Matlab函数四、Matlab文件操作初步求微分方程的解一、问题背景和实验目的二、求解微分方程的若干计算方法三、求解微分方程的Matlab函数四、Matlab文件操作初步一、问题背景和实验目的一、问题背景和实验目的实验四、求微分方程的解实验四、求微分方程的解问题背景和实验目的自牛顿发明微积分以来微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程绝大多数是微分方程。由于实际应用的需要人们必须求解微分方程。然而能够求得解析解的微分方程十分有限绝大多数微分方程需要利用数值方法来近似求解。本实验主要研究如何用Matlab来计算微分方程(组)的数值解并重点介绍一个求解微分方程的基本数值解法--Euler折线法。微分方程模型举例微分方程模型举例导弹动力学弹道方程组(垂直平面)在弹道设计中求解动力学弹道的目的是为了得到xyv三个参数以便对射程、导引方法及燃料添置等方案进行选择。微分方程模型举例微分方程模型举例金融工程期权定价模型年Black和Scholes按下面的微分方程决定欧式期权的无套利价格在基于股票的衍生证券市场上欧式买入期权的行使办法是:在到期日T当股票价格ST>X(行使价格)时则按X买入股票否则不行使期权。欧式卖出期权的行使办法是:在到期日T当股票价格ST<X(行使价格)时则按X卖出股票否则不行使期权。其中r为无风险利率S为股票价格X为期权的行使价格为基于S的期权价格。二、求解微分方程的若干计算方法二、求解微分方程的若干计算方法初值问题的Euler折线法初值问题的Euler折线法考虑一维经典初值问题基本思想:用差商代替微商根据Talyor公式y(x)在点xk处有==>初值问题的Euler折线法(续)初值问题的Euler折线法(续)具体步骤:分割求解区间差商代替微商解代数方程。给定等距剖分:步长为分割求解区间差商代替微商==>得方程组:Euler折线法举例Euler折线法举例例:用Euler法解初值问题解:取步长h=()n=n得差分方程当h=即n=时Matlab源程序见附录。附录:Euler折线法源程序附录:Euler折线法源程序clearf=sym('y*xy^')a=b=h=n=(ba)hn=(ba)hx=y=szj=x,yfori=:ni=:ny=yh*subs(f,{'x','y'},{x,y})x=xhszj=szjx,yendszjplot(szj(:,),szj(:,),'or')Euler折线法举例(续)Euler折线法举例(续)可以求得解析解为:解析解近似解四阶RungeKutta方法四阶RungeKutta方法由图可知近似解的误差较大为了减小误差可采用以下方法:让步长h取得更小一些改用具有较高精度的数值方法:RungeKutta(龙格-库塔)方法在实际应用中用得最多的是四阶RK方法(教材第页)附录:四阶RK方法源程序附录:四阶RK方法源程序clearf=sym('y*xy^')a=b=h=n=(ba)hn=(ba)hx=y=szj=x,yfori=:ni=:nl=subs(f,{'x','y'},{x,y})l=subs(f,{'x','y'},{xh,yl*h})l=subs(f,{'x','y'},{xh,yl*h})l=subs(f,{'x','y'},{xh,yl*h})y=yh*(l*l*ll)x=xhszj=szjx,yendplot(szj(:,),szj(:,),'dg')四阶RungeKutta方法四阶RungeKutta方法Euler法与RK法误差比较Euler法与RK法误差比较三、求解微分方程的Matlab函数三、求解微分方程的Matlab函数用Matlab函数解初值问题用Matlab函数解初值问题用Maltab函数解初值问题求解析解:dsolve求数值解:ode、ode、ode、odetodes、odes、odetbdsolve使用格式y=dsolve('eq','eq',,'cond','cond',,'v')求微分方程的解析解其中y为输出eq、eq、为微分方程cond、cond、为初值条件v为自变量。例:求微分方程的通解并加以验证。y=dsolve('Dy*x*y=x*exp(x^)','x')dsolve的使用dsolve的使用几点说明y=dsolve('eq','eq',,'cond','cond',,'v')例:y=dsolve('Dy*x*y=x*exp(x^)','x')、微分方程中用D表示对自变量的导数如:、如果省略初值条件则表示求通解、如果省略自变量则默认自变量为tdsolve('Dy=*x','x')%dydx=xdsolve('Dy=*x')%dydt=x、若找不到解析解则返回其积分形式。dsolve举例(例)dsolve举例(例)例:求微分方程在初值条件下的特解并画出解函数的图形。y=dsolve('x*Dyyexp(x)=','y()=*exp()','x')ezplot(y)dsolve举例(例)dsolve举例(例)例:求微分方程组在初值条件下的特解并画出解函数的图形。x,y=dsolve('Dx*xy=exp(t)','Dyx*y=','x()=','y()=','t')ezplot(x,y,,)注:解微分方程组时如果所给的输出个数与方程个数相同则方程组的解按词典顺序输出如果只给一个输出则输出的是一个包含解的结构(structure)类型的数据。dsolve举例(例)dsolve举例(例)例:x,y=dsolve('Dx*x=','Dy*y=','x()=','y()=','t')r=dsolve('Dx*x=','Dy*y=','x()=','y()=','t')这里返回的r是一个结构类型的数据rx查看解函数x(t)ry查看解函数y(t)只有很少一部分微分方程(组)能求出解析解。大部分微分方程(组)只能利用数值方法求数值解。dsolve的输出个数只能为一个或与方程个数相等。初值问题数值解初值问题数值解用Maltab函数求初值问题的数值解T,Y=solver(odefun,tspan,y)求微分方程的数值解。其中y为初值条件tspan为求解区间系统在数值求解时自动对求解区间进行分割T(向量)中返回的是分割点的值(自变量)Y(向量)中返回的是解函数在这些分割点上的函数值。solver为Matlab的ODE求解器(可以是ode、ode、ode、odes、odes、odet、odetb)因为没有一种算法可以有效地解决所有的ODE问题为此MATLAB提供了多种求解器solver对于不同的ODE问题采用不同的solver。ODE求解器的特点ODE求解器的特点不同ODE求解器solver的特点参数说明参数说明T,Y=solver(odefun,tspan,y)odefun为显式常微分方程可以用命令inline定义或在函数文件中定义然后通过函数句柄调用。fun=inline('*y*x^*x','x','y')x,y=ode(fun,,,)注:也可以在tspan中指定对求解区间的分割如:x,y=ode(fun,::,)ODE数值求解举例ODE数值求解举例如果需求解的问题是高阶常微分方程则需将其化为一阶常微分方程组此时需用函数文件来定义该常微分方程组。ODE数值求解举例ODE数值求解举例先编写函数文件verderpolmfunctionxprime=verderpol(t,x)globalmuxprime=x(),mu*(x()^)*x()x()再编写命令文件vdplm在命令窗口直接运行该文件。clearglobalmuy=t,x=ode('verderpol',,,y)t,x=ode(verderpol,,,y)plot(t,x(:,),'r',t,x(:,),'b')关于ODE函数定义的进一步说明关于ODE函数定义的进一步说明用Matlab函数ode、ode等计算微分方程(组)的数值解时如何定义odefun?Matlab中如何定义函数(ODE)如何定义函数(ODE)如何定义函数(ODE)T,Y=ode(odefun,tspan,y)T,Y=ode(odefun,tspan,y)ode、ode等函数用于求解显式常微分方程当是向量函数时所对应的上面的方程即为微分方程组odefun举例说明举例说明fun=inline('*y*x^*x','x','y')x,y=ode(fun,,,)解法一:使用inline定义微分方程odefunodefun为方程右端项f(t,y)可以用inline定义(只适合于单个方程的情形)通过函数文件定义然后用函数句柄调用(适合所有情形)举例说明(单个方程)举例说明(单个方程)functiondy=myfun(x,y)dy=*y*x^*x解法二:通过函数文件定义微分方程odefun、先编写函数文件myfunmclearx,y=ode(myfun,,,)、编写主文件mainm或直接在Matlab命令窗口输入上面的语句。举例说明(方程组)举例说明(方程组)解:此时只能通过函数文件定义微分方程odefun例:求,,的数值解。functiondy=myfun(t,y)dy=zeros(,)dymustbeacolumnvector!dy()=y()*y()dy()=y()*y()dy()=*y()*y()、先编写函数文件myfunmclearT,Y=ode(myfun,,,,,)、编写主文件mainmdy=y()*y()y()*y()*y()*y()思考思考functiondy=myfun(t,y)dy=zeros(,)dy()=y()*y()dy()=y()*y()dy()=*y()*y()、函数文件myfunm能不能写成下面形式?functiondy=myfun(t,x,y,z)dy=zeros(,)dy()=y*zdy()=x*zdy()=*x*yX说明():说明():odefun变量属性必须一一对应!functiondy=myfun(t,y)如果是常微分方程组y就是列向量!dy必须是列向量长度为方程组的个数通常与y的长度相同!函数中的输入参数和输出参数是形参名字可以任意取但必须满足上述条件。即输入参数有两个第一个表示自变量第二个是由应变量组成的列向量输出参数必须是列向量。说明():说明():T,Y=ode(odefun,tspan,y)T,Y=ode(odefun,tspan,y)odefun?ode、ode能解什么样的ODE?例例functiondy=myfun(t,y)dy=zeros(,)dy()=y()tdy()=t例:解初值问题:functionyprime=myfun(t,y)yprime=y()tt–functionyprime=myfun(x,y)yprime=y()xx–clearT,Y=ode(myfun,,,,)、主文件mainm、函数文件myfunm高阶常微分方程高阶常微分方程高阶常微分方程例:VerderPol初值问题一阶常微分方程组变量替换化为参数怎么处理?用全局变量传递参数传递参数传递、函数文件verderpolmfunctionxprime=verderpol(t,x)globalmuxprime=x()mu*(x()^)*x()x()、主文件vdpmclearglobalmuy=mu=t,x=ode('verderpol',,,y)t,x=ode(verderpol,,,y)plot(t,x(:,),'r',t,x(:,),'b')四、Matlab文件操作初步四、Matlab文件操作初步用Matlab函数dsolve求出微分方程(组)的解析解后如何将这个解析解打印到文件中?如何将符号对象输出到文件中如何将符号对象输出到文件中用Matlab函数dsolve求出微分方程(组)的解析解后如何将这个解析解打印到文件中?例:求微分方程的通解并将其输出到文件solutiontxt中。y=dsolve('Dy*x*y=x*exp(x^)','x')whosy此时y是一个符号对象fprintf无法打印一个符号对象需要先用char函数将其转化成一个字符串后再按字符串格式输出。fprintf('Thesolutionis:y=sn',char(y))char一个完整的程序一个完整的程序cleary=dsolve('Dy*x*y=x*exp(x^)','x')下面的语句是在屏幕上输出结果fprintf('Thesolutionis:y=sn',char(y))下面是将结果输出到指定的文件fid=fopen('solutiontxt','wt')fprintf(fid,'TheSolutionis:y=sn',char(y))fclose(fid)微分方程组举例微分方程组举例例:求微分方程组在初值条件下的特解并将其追加到文件solutiontxt中。clearx,y=dsolve('Dx*xy=exp(t)','Dyx*y=','x()=','y()=','t')fid=fopen('solutiontxt','at')fprintf(fid,'nSolution:nx=sny=sn',char(x),char(y))fclose(fid)只给一个输出时的情形只给一个输出时的情形clearsol=dsolve('Dx*xy=exp(t)','Dyx*y=','x()=','y()=','t')fid=fopen('solutiontxt','at')fprintf(fid,'nSolution:nx=sny=sn',char(solx),char(soly))fclose(fid)作业作业教材P:练习~要求:将练习、的结果(数值解)分别写入文件solextxt和solextxt

用户评价(0)

关闭

新课改视野下建构高中语文教学实验成果报告(32KB)

抱歉,积分不足下载失败,请稍后再试!

提示

试读已结束,如需要继续阅读或者下载,敬请购买!

评分:

/44

VIP

在线
客服

免费
邮箱

爱问共享资料服务号

扫描关注领取更多福利