下载

1下载券

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

上传资料

关闭

关闭

关闭

封号提示

内容

首页 MATLAB解微分方程

MATLAB解微分方程.ppt

MATLAB解微分方程

风凌万里
2011-08-23 0人阅读 举报 0 0 暂无简介

简介:本文档为《MATLAB解微分方程ppt》,可适用于高等教育领域

MATLABODE初值问题的数值解PDE问题的数值解问题提出倒葫芦形状容器壁上的刻度问题对于如图所示圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式其中直径D为常数对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢下表是经过测量得到部分容器高度与直径的关系xo根据上表的数据,可以拟合出倒葫芦形状容器的图,建立如图所示的坐标轴,问题为如何根据任意高度x标出容器体积V的刻度,由微元思想分析可知HD只要求解上述方程,就可求出体积V与高度x之间的函数关系,从而可标出容器壁上容积的刻度,但问题是函数D(x)无解析表达式,无法求出其解析解因此,得到如下微分方程初值问题包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。微分方程分类常微分方程:()()式称为初值问题在实际应用中还经常需要求解常微分方程组:()式称为边值问题。但能求解析解的常微分方程是有限的大多数的常微分方程是给不出解析解的这个一阶微分方程就不能用初等函数及其积分来表达它的解。例例事实上从实际问题当中抽象出来的微分方程通常主要依靠数值解法来解决。微分方程数值方法的基本思想对常微分方程初值问题的数值解法就是要算出精确解y(x)在区间a,b上的一系列离散节点处的函数值相邻两个节点的间距称为步长步长可以相等也可以不等。假定h为定数称为定步长这时节点可表示为数值解法需要把连续性的问题加以离散化从而求出离散节点的数值解。离散化对常微分方程数值解法的基本出发点就是离散化。其数值解法的基本特点:采用“步进式”即求解过程顺着节点排列的次序一步一步地向前推进描述这类算法要求给出用已知信息计算的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法对初值问题数值解和精确解数值解和精确解用数值方法求解初值问题不是求出它的解析解或其近似解析式而是给出它的解在某些离散节点上的近似值。用y(x)表示问题的准确解y(x),y(x),y(xN)表示解y(x)在节点x,x,…,xN处的准确值y,y,…,yN表示数值解即问题的解y(x)在相应节点处的近似值。单步法和多步法单步法和多步法单步法:在计算yi时只利用yi多步法:在计算yi时不仅利用yi,还要利用yi−,yi−,…,k步法:在计算yi时要用到yi,yi−,…,yi−k显式格式可写成:yk=ykhΦf(xkykh)隐式格式:yk=ykhΦf(xkyk,ykh)它每步求解yk需要解一个隐式方程。欧拉(Euler)方法在x=x处用差商代替导数:由得同理在x=xn处用差商代替导数:由得若记则上式可记为此即为求解初值问题的Euler方法又称显式Euler方法。例:用Euler方法求解常微分方程初值问题并将数值解和该问题的解析解比较。解:Euler方法的具体格式:h=y()=x=:h:forn=:xn=x(n)yn=y(n)y(n)=ynh*(ynxn*yn*yn)endx=:h:y=x(x^)plot(x,y,x,y,x,y,'o')程序实现xny(xn)ynyny(xn)h=,xn=nh,(n=,,…,),f(x,y)=yx–y计算中取f(,)=计算结果如下:xny(xn)ynyny(xn)由表中数据可以看到微分方程初值问题的数值解和解析解的误差一般在小数点后第二位或第三位小数上说明Euler方法的精度是比较差的。二阶RungeKutta方法其中c,c,,待定。上式的局部截断误差为:即c=a==(a)方程组解不唯一可令c=a则满足上述条件的公式都为阶RK公式。例蛇形曲线的初值问题令f(x,y)=yx–y,取f(,)=,h=,xn=hn,(n=,,…,)阶龙格库塔公式计算格式:k=ynxn–yn,k=(ynhk)(xnh)–(ynhk)yn=ynhkkx=y=h=x=:h:k=k=(yh*k)x()*(yh*k)^y()=y*h*(kk)forn=:k=y(n)x(n)*y(n)^k=(y(n)h*k)x(n)*(y(n)h*k)^y(n)=y(n)*h*(kk)endy=x(x^)plot(x,y,'o',x,y)常用的一个公式为四阶RungeKutta方法functionydot=harmonic(t,y)ydot=y()y()y=inline(‘*y’,’t’,’y’)SystemofEquationsfunctionydot=twobody(t,y)r=sqrt(y()^y()^)ydot=y()y()y()r^y()r^TwoBodyProblemLinearizedDifferentialEquationsJ的特征值是解增长解衰减解振荡MATLAB求常微分方程数值解的函数MATLAB求常微分方程数值解的函数基于龙格-库塔法MATLAB求常微分方程数值解的函数一般调用格式为:t,y=ode('fname',tspan,y)t,y=ode('fname',tspan,y)其中fname是定义f(t,y)的函数文件名该函数文件必须返回一个列向量。tspan形式为t,tf,表示求解区间。y是初始状态列向量。t和y分别给出时间向量和相应的状态向量。ode:Bogacki,Shampine()和Shampine(),””表示用两算法:一个阶一个阶Bogacki,PandLFShampine,"A()pairofRungeKuttaformulas,"ApplMathLetters,Vol,,ppBSalgorithmF=inline('y()y()','t','y')ode(F,*pi,)opts=odeset(‘reltol’,e,’abstol’,e,’outpcn’)Examplesode(twobody,*pi,)Examplesy=ode(twobody,*pi,y)y=t,y=ode(twobody,*pi,y)plot(y(:,),y(:,))axisequaly=t,y=ode(twobody,*pi,y)plot(y(:,))plot(y(:,))Aproblemisstiffifthesolutionbeingsoughtisvaryingslowly,buttherearenearbysolutionsthatveryrapidly,sothenumericalmethodmusttakesmallstepstoobtainsatisfactoryresultsAmodelofflamepropagation(火焰燃烧):y是球的半径y^和y^与球的表面积和体积有关想一下点燃一根火柴光球迅速增长到达一个关键的大小然后维持它的大小(由于进入球内氧气和消耗的氧气平衡)StiffProblem(刚性问题)eta=symyF=inline(‘y^y^’,’t’,’y’)ode(F,eta,eta)eta=symyF=inline(‘y^y^’,’t’,’y’)ode(F,eta,eta)eta=odes(inline('y^y^','t','y'),eta,eta)例蛇形曲线的常微分方程初值问题MATLAB数值求解命令F=inline('(x^)*y^')ode(F,,,)输出结果为图形T,y=ode(f,,,)将得到自变量和函数的离散数据T=y=例洛伦兹模型由如下常微分方程组描述取===。初值:x()=y()=z()=。利用MATLAB求解常微分方程数值解命令计算出t∈内三个未知函数的数据值并绘出相空间在YX平面的投影曲线气象学家Lorenz提出一篇论文名叫「一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风?」论述某系统如果初期条件差一点点结果会很不稳定他把这种现象戏称做「蝴蝶效应」。Lorenz为何要写这篇论文呢?这故事发生在年的某个冬天他如往常一般在办公室操作气象电脑。平时他只需要将温度、湿度、压力等气象数据输入电脑就会依据三个内建的微分方程式计算出下一刻可能的气象数据因此模拟出气象变化图。这一天Lorenz想更进一步了解某段纪录的後续变化他把某时刻的气象数据重新输入电脑让电脑计算出更多的後续结果。当时电脑处理数据资料的数度不快在结果出来之前足够他喝杯咖啡并和友人闲聊一阵。在一小时後结果出来了不过令他目瞪口呆。结果和原资讯两相比较初期数据还差不多越到後期数据差异就越大了就像是不同的两笔资讯。而问题并不出在电脑问题是他输入的数据差了而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。天气预报的准确性:http:wwwnanicomtwseniorseniorteachstfastnewsstnfastnewsstnfastnewsstnfastnewshtmLorenz现象的数学:http:wwwsciscapeorgnewsdetailphpnewsid=分形艺术电子版:http:wwwphilpkueducnpersonalhuajiefractalarthtmlnoticehtm混沌理论:http:wwwlamostorg~yzhaodocchaoshtml记向量yyy=xyz创建MATLAB函数文件如下functionz=flo(t,y)z(,:)=*y()y()*y()z(,:)=*(y()y())z(,:)=y()*y()*y()y()用MATLAB命令求解并绘出YX平面的投影图y=x,y=ode(flo,,,y)plot(y(:,),y(:,))plot(y(:,),y(:,))plot(y(:,),y(:,),y(:,))y=plot(y(:,i))非刚性系统:ode(RungeKutta)ode(RungeKuatta)ode(AdamsBashforthMoultonPECE)多步方法刚性系统:odes(Gear方法),多步方法odes(二阶modifiedRosenbroackformula)单步odet(trapezoidalrule),solveDAEsodetb(TRBDF)lowordermethodMatlab'sODESolversLaplacian算子:Poisson方程(elliptic):Laplacian算子的特征值问题:Heatequation(parabolic):Waveequation(hyperbolic):PDEModel五点离散Poisson方程离散:特征值问题:FineteDifferenceMethods热方程:波动方程:Matrixs椭圆方程:特征值方程:热方程:波动方程:离散方程例用古典显式格式求解抛物型方程初始条件为边界条件为取步长⊿x=h=,⊿t=k=。解r=kh==,古典显式格式为,,,,,┇,,┇,functionu=gudian(f,a,b,c,c,m,n)输入初值和Uh=a(m)k=b(n)r=kh^U=zeros(n,m)赋边界条件U(:n,)=cU(:n,m)=cfunctiony=fg(x)y=*(xx^)赋初始条件U(,:m)=fg(:h:h*(m))计算内点上u的数值解Ufori=:nforj=:(m)U(i,j)=(*r)*U(i,j)r*(U(i,j)U(i,j))endendgudianlm步长h=,k=,r=kh=a=b=c=c=m=n=U=gudian('fg',a,b,c,c,m,n)x=::ay=::bX,Y=meshgrid(x,y)surf(X,Y,U)输入U后再画图PDETOOL有限元方法

用户评价(0)

关闭

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

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

提示

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

文档小程序码

使用微信“扫一扫”扫码寻找文档

1

打开微信

2

扫描小程序码

3

发布寻找信息

4

等待寻找结果

我知道了
评分:

/65

MATLAB解微分方程

VIP

在线
客服

免费
邮箱

爱问共享资料服务号

扫描关注领取更多福利