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

上传资料

关闭

关闭

关闭

封号提示

内容

首页 Matlab解微分方程[整理版]

Matlab解微分方程[整理版].doc

Matlab解微分方程[整理版]

himikokiki
2017-10-22 0人阅读 举报 0 0 暂无简介

简介:本文档为《Matlab解微分方程[整理版]doc》,可适用于职业岗位领域

Matlab解微分方程整理版第十六章,,,,,,,,,,偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程称为偏微分方程(partial,,,,,differential,,,,,equation)简称PDE。偏微分方程问题其求解是十分困难的。除少数特殊情况外绝大多数情况均难以求出精确解。因此近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。,,,,,,,,,,几类偏微分方程的定解问题一个偏微分方程的表示通常如下:,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,ABCfxy,,,,,,,(,,,,)xxxyyyxy,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()式中是常数称为拟线性(quasilinear)数。通常存在种拟线性方程:ABC,,BAC,,,,,,,,,,,,,,,,,,,,,,双曲型(hyperbolic)方程:BAC,,,,,,,,,,,,,,,,,,,,,,抛物线型(parabolic)方程:BAC,,,,,,,,,,,,,,,,,,,,,,椭圆型(ellliptic)方程:。,,,,,,,,,,,,,,,双曲型方程最简单形式为一阶双曲型方程:,,uu,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,tx,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()物理中常见的一维振动与波动问题可用二阶波动方程:,,uu,,,,,,,,,,,,,,'g',,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,and,,,,,b,,,,,right,,,,,endpoints,,,,,of,,,,,,a,,,,,and,,,,,,b,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,c,,,,,the,,,,,constant,,,,,in,,,,,the,,,,,wave,,,,,equation,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,n,,,,,and,,,,,m,,,,,number,,,,,of,,,,,grid,,,,,points,,,,,over,,,,,,a,,,,,and,,,,,,bOutput,,,,,,,,,,U,,,,,solution,,,,,matrix,,,,,analogous,,,,,to,,,,,Table,,,,,Initialize,,,,,parameters,,,,,and,,,,,U,,,,,h,,,,,=,,,,,a(n)k,,,,,=,,,,,b(m)r,,,,,=,,,,,c*khr=r^r=r^s,,,,,=,,,,,,,,,,,,,,,r^s,,,,,=,,,,,,,,,,,,,,,*r^U,,,,,=,,,,,zeros(n,m),,,,,Comput,,,,,first,,,,,and,,,,,second,,,,,rows,,,,,for,,,,,i=:n,,,,,,,,,,,,,,,U(i,)=feval(f,h*(i)),,,,,,,,,,,,,,,U(i,)=s*feval(f,h*(i))k*feval(g,h*(i)),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,r*(feval(f,h*i)feval(f,h*(i)))end,,,,,,,,,,,,,,,,,,,,Compute,,,,,remaining,,,,,rows,,,,,of,,,,,U,,,,,for,,,,,j=:m,,,,,,,,,,,for,,,,,i=:(n),,,,,,,,,,,,,,,,,,,,,,,,,,U(i,j),,,,,=,,,,,s*U(i,j)r*(U(i,j)U(i,j))U(i,j),,,,,,,,,,endend,,,,,U=U'*******************************************************************,,,,,,,,,,抛物型方程的差分解法以一维热传导方程:,,uu,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,aa(),,tx,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()为基本模型讨论适用于抛物型方程定解问题的几种差分格式。,,,,,,,,,,差分格式的建立首先对平面进行网格剖分。分别取为x方向与方向的步长用两族平行直线h,,xt,tttjj,,,,(,,)?xxkhk,,,,,(,,,)?将平面剖分成矩形网格节点为:xt,jk(,)(,,,,,,)xtkj,,,,??。为简便记:kj(,)(,)kjxt,ukjuxt(,)(,),ggt,(),,,()x,,,()x,,,,,,,,,,,,,,,,,,,,kjkjjjkkkk。,,,()t,,,()tjjjj,,,,,,,,,,微分方程的差分近似,u在网格内点处对分别采用向前、向后及中心差商公式:(,)kj,t,,uukjukj(,)(,),,O(),t,(,)kj,,,uukjukj(,)(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,O(),,t,(,)kj,,,uukjukj(,)(,),O(),,t,(,)kj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()一维热传导方程可分别表示为:ukjukjukjukjukj(,)(,)(,)(,)(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,aOh()h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()ukjukjukjukjukj(,)(,)(,)(,)(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,aOh()h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()ukjukjukjukjukj(,)(,)(,)(,)(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,aOh(),,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()由此得到一维热传导方程的不同差分近似:uuuuu,,kjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()uuuuu,,kjkjkjkjkj,,,,,,,,,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()uuuuu,,kjkjkjkjkj,,,,,,,,,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()上述差分方程所用到的节点各不相同。其截断误差分别为和Oh(),Oh(),。因此它们都与一维热传导方程相容。Oh(),ukjukjukjukjukj(,)(,)(,)(,)(,),,,,如果将式:中的用,u,,aOh()kj,,h代替则可得到又一种差分近似:uu()kjkj,,,uuuuuu,,,kjkjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()差分方程用到四个节点。由Taylor公式容易得出:uuuO,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()()kjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(),,,故其的截断误差为OhO()。因而不是对任意的此差分方程都能逼近热h,,,,,,h,,,,uu传导方程:仅当时才成立。,,,aa,(),,oh(),,tx综上可知用不同的差商公式可以得到微分方程的不同的差分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳定性。前面三个方程都具有相容性而此方程则要在一定条件下才具有相容性。,,,,,,,,,,初、边值条件的处理为用差分方法求解定解问题初值问题:,,,uu,,,,,,,,atx,,,,tx,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uxxx(,)(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()初边值问题:,,,uuatTxl,,,,,,,,,,tx,,,uxxxl(,)(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,utgtultgttT(,)(),(,)(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()还需对定解条件进行离散化。对初始条件及第一类边界条件可直接得到:,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uux,,(,),(,,,,,)kkn,,,??或kkk,uutg,,(,),jjj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,,,)jm,,?uultg,,(,)njjj,lTnm,,,式中:。,h对第二、三类边界条件:,u,,,()()tugtx,,,,x,,tT,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,u()(),tugt,xl,,,,x,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()需用差分近似。下面介绍两种较简单的处理方法。,u,u在左边界处用向前差商近似偏导数在右边界处用向后差商近似()x,()xl,,x,x即:,,uujuj(,)(,),Oh(),xh(,)j,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uunjunj(,)(,),Oh(),xh(,)nj(,,,)jm,?,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()则得边界条件的差分近似为:uu,,,,jj,,,ug,jjj,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uu,,njnj,,,,ug,,jnjj,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,,,)jm,?,,,,,,,,,()其截断误差为。Oh(),u(),,,,,用中心差商近似,即:,x,,,uujuj(,)(,),Oh(),xhj(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uunjunj(,)(,),Oh(),xhnj(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()(,,,)jm,?则得边界条件的差分近似为:uu,,,,jj,,,,ug,jjj,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,,,)jm,?,uu,njnj,,,,,ug,,jnjj,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()其截断误差为。误差的阶数提高了但出现定解区域外的节点和这就(,),j(,)njOh()uu需要将解拓展到定解区域外。可以通过用内节点上的值插值求出和也可以假定u,,jnj,u热传导方程在边界上也成立将差分方程扩展到边界节点上由此消去和。un,j,,j,,,,,几种常用的差分格式以热传导方程的初边值问题:,,,uu,,,,,,atTxl,,,,tx,,uxxxl(,)(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,utgtultgttT(,)(),(,)(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()为例给出几种常用的差分格式。(),,,,,古典显式格式a,令则:r,huuuuu,,kjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()可改写成:,,,,,。将其与初始条件及第一类边界条件:urururu,,()kjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uux,,(,),(,,,,,)kkn,,,??或kkk,uutg,,(,),jjj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,,,)jm,,?uultg,,(,)njjj,结合我们得到求解此问题的一种差分格式:,urururuknjm,,,,,,()(,,,,,,,??kjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,ukn,,(,,,)?,kk,,ugugjm,,,,(,,,)?,,jjnjj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()(u,),由于第层上节点处的值已知由此即可算出在第一层上节点处uu()j,()j,k,kuu的近似值。重复使用此式可以逐层计算出所有的因此此差分格式称为典显式格式。k,kj,又因式中只出现相邻两个时间层的节点故此式是二层显式格式。(),,,,,古典隐式格式将式:uuuuu,,kjkjkjkjkj,,,,,,,,,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()整理并与初始条件及第一类边界条件式联立得差分格式如下:,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uuruuuknjm,,,,,,()(,,,,,,,)??kjkjkjkjkj,,,,,,,,,,,,,,,,,(),ukn,,(,,,)?,kk,,ugugjm,,,,(,,,)?,,jjnjj,a,其中:。虽然第层上的值仍为已知但不能由上式直接计算以上各层节点上的值r,uukj,h必须通过解下列线性方程组:,rr,,uurg,,rrr,,jjj,,,,,,uu,,,,jj,,,,,,,,,,??,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uunjnj,,,,,,,,,,,,,,uurg,,rrrnjnjj,,,,,,,,,rr,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()uu式中。,,,,,才能由计算故此差分格式称为古典隐式格式。此方程(,,,)jm,,?kj,kj,组是三对角方程组且系数矩阵严格对角占优故解存在唯一。(),,,,,Richardson格式Richardson格式是将式:uuuuu,,kjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,h整理后与初始条件及第一类边界条件式联立。其计算公式为:,uuruuuknjm,,,,,,()(,,,,,,,??kjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,ukn,,(,,,)?,kk,,ugugjm,,,,(,,,)?,,jjnjj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,)j,,j,j这种差分格式中所涉及的节点出现在三层上故为三层显式格式。Richardson格式是一种完全不稳定的差分格式因此它在实际计算中是不能采用的。(),,,,,杜福特弗兰克尔,,,,,(DoFortFrankel),,,,,格式DoFortFrankel格式也是三层显式格式它是由式:uuuuuu,,,kjkjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,h与初始条件及第一类边界条件式结合得到的。具体形式如下:,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,rr,,uuuuknjm,,,,,()(,,,,,,,)??kjkjkjkj,,,,,,,rr,,,,,,,,,,,(),ukn,,(,,,)?,kk,,ugugjm,,,,(,,,)?,,jjnjj,,用这种格式求解时除了第层上的值由初值条件得到必须先用二层格式求出第层上uk,的值然后再按上式逐层计算。uujm(,,,),?k,kj,(),,,,,六点隐式格式对二阶中心差商公式:,,,uukjukjukj(,)(,)(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,xhkj(,),u,u(,)kj如果用在与处的二阶中心差商的平均值近似在处的值即:(,)kj(,)kj,x,xuuuuuu,,,,ukjkjkjkjkjkj,,,,,,,,,Oh(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,xhkj(,),u(k,j)同时在点处的值也用中心差商:,t,,,uukjukj(,)(,),,O(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,(,)kj近似即:uu,,ukjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,x(,)kj这样又得到热传导方程的一种差分近似:uu,akjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,),,,,uuuuuukjkjkjkjkjkj,,,,,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()其截断误差为将上式与初始条件及第一类边界条件式联立并整理得差分格式:O(,h)rr,uuuuuuuu,,,()()kjkjkjkjkjkjkjkj,,,,,,,,,,,,,(,,,,,,,)knjm,,,,??,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,ukn,,(,,,)?kk,,ugugjm,,,,(,,,)?,,,jjnjj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()此格式涉及到六个节点它又是隐式格式故称为六点隐式格式。与古典隐式格式类似用六点格式由第层的值计算第层的值时需求解三对角方程组:uujjkj,kj,,rr,,u,,rrr,j,,,,u,,,j,,,,,,?,,,,,,unj,,,,,,,,u,,rrrnj,,,,,,,rrrruuuuu,(),,,,,jjjjj,,,,r,,uuuu,(),,,,jjjj,,,,,?,,,,ruuuu,()njnjnjnj,,,,,,,,,,,,rr,,u,uuuu()n,,,,,,jnjnjnjnj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()式中(,,,)jm,,?。此方程组的系数矩阵严格对角占优故仍可用追赶法求解。例,,,,,,,,,,用古典显式格式求初边值问题。,,,uu,,,,,,,tx,,,tx,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uxxx(,),,,,,ututt(,),(,),,,,,,,的数值解取。h,,,,a,解:这里:。ar,,,,gtgt(),(),,,()xx,h由格式:,urururuknjm,,,,,,()(,,,,,,,??kjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,ukn,,(,,,)?,kk,,ugugjm,,,,(,,,)?,,jjnjj,可得到:uuukj,,,()(,,,,,)?,kjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uxkk,,,(,,,),kk,,uuj,,,,(,,,)?,,jj,将初值u代入上式即可算出:k,uuu,,,()(),,,uuu,,,()(),,,u,u,将边界条件及上述结果代入又可求得:,,uuuu,,,,,,,,,,,如此逐层计算得全部节点上的数值解为:uuuu,,,,,,,,,,,,u,,??,,,,,,,,,,,前向差分法求解热传导方程的MATLAB程序,,xa,,tbutcuatc(,),(,),,,,,,,,,,,,,,,,,,,,,,设uxfx(,)(),其中而且其中Rxtxatb,,,,,{(,):,}uxtcuxt(,)(,),求解区间内的近似解。txx*****************************************************function,,,,,U=forwdif(f,c,c,a,b,c,n,m),,,,,Input,,,,,,,,,,f=u(x,),,,,,as,,,,,a,,,,,string,,,,,'f',,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,c=u(,t),,,,,and,,,,,c=u(a,t),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,and,,,,,b,,,,,right,,,,,endpoints,,,,,of,,,,,,a,,,,,and,,,,,,b,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,c,,,,,the,,,,,constant,,,,,in,,,,,the,,,,,heat,,,,,equation,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,n,,,,,and,,,,,m,,,,,number,,,,,of,,,,,grid,,,,,points,,,,,over,,,,,,a,,,,,and,,,,,,bOutput,,,,,,,,,,U,,,,,solution,,,,,matrix,,,,,analogous,,,,,to,,,,,Table,,,,,Initialize,,,,,parameters,,,,,and,,,,,U,,,,,h=a(n)k=b(m)r=c^*kh^s=*rU=zeros(n,m),,,,,Boundary,,,,,conditions,,,,,U(,:m)=cU(n,:m)=c,,,,,Generate,,,,,first,,,,,rowU(:n,)=feval(f,h:h:(n)*h)',,,,,Generate,,,,,remaining,,,,,rows,,,,,of,,,,,U,,,,,for,,,,,j=:m,,,,,,,,,,,,,,,for,,,,,i=:n,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,U(i,j)=s*U(i,j)r*(U(i,j)U(i,j)),,,,,,,,,,,,,,,endend,,,,,U=U'*****************************************************,,,,,,,,,,CrankNicholson求解热传导方程的MATLAB程序,,xa,,tb,,,,,,,,,,,,,,,,,,,,设其中而且其中utcuatc(,),(,),,uxfx(,)(),求解区间内的近似解。Rxtxatb,,,,,{(,):,}uxtcuxt(,)(,),txx*****************************************************function,,,,,U=crnich(f,c,c,a,b,c,n,m),,,,,Input,,,,,,,,,,f=u(x,),,,,,as,,,,,a,,,,,string,,,,,'f',,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,c=u(,t),,,,,and,,,,,c=u(a,t),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,and,,,,,b,,,,,right,,,,,endpoints,,,,,of,,,,,,a,,,,,and,,,,,,b,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,c,,,,,the,,,,,constant,,,,,in,,,,,the,,,,,heat,,,,,equation,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,n,,,,,and,,,,,m,,,,,number,,,,,of,,,,,grid,,,,,points,,,,,over,,,,,,a,,,,,and,,,,,,bOutput,,,,,,,,,,U,,,,,solution,,,,,matrix,,,,,analogous,,,,,to,,,,,Table,,,,,Initialize,,,,,parameters,,,,,and,,,,,U,,,,,h=a(n)k=b(m)r=c^*kh^s=rs=rU=zeros(n,m),,,,,Boundary,,,,,conditions,,,,,U(,:m)=cU(n,:m)=c,,,,,Generate,,,,,first,,,,,row,,,,,U(:n,)=feval(f,h:h:(n)*h)',,,,,Form,,,,,the,,,,,diagonal,,,,,and,,,,,offdiagonal,,,,,elemnts,,,,,of,,,,,A,,,,,and,,,,,the,,,,,constant,,,,,vector,,,,,B,,,,,and,,,,,solve,,,,,tridiagonal,,,,,system,,,,,AX=B,,,,,Vd(,:n)=s*ones(,n)Vd()=Vd(n)=Va=ones(,n)Va(n)=Vc=ones(,n)Vc()=Vb()=cVb(n)=cfor,,,,,j=:m,,,,,,,,,,,,,,,for,,,,,i=:n,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,Vb(i)=U(i,j)U(i,j)s*U(i,j),,,,,,,,,,,,,,,end,,,,,,,,,,,,,,,X=trisys(Va,Vd,Vc,Vb),,,,,,,,,,,,,,,U(:n,j)=X'end,,,,,U=U'*****************************************************,,,,,,,,,,椭圆型方程第一边值问题的差分解法本节以Poisson方程为基本模型讨论第一边值问题的差分方法。,,,,,,,,,,差分格式的建立考虑Poisson方程的第一边值问题:,,,uu(,)(,),,,fxyxy,,,xy,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,)(,)uxyxy,,,,,xy,,(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()图取分别为方向和方向的步长如图所示以两族平行线:yxh,,yyjkj,,,,,,(,,,,)?xxkh,,jk将定解区域剖分成矩形网格。节点的全体记为:Rxyxkhyjkj,,,{(,),,,},为整数,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,。kjkjR:,定解区域内部的节点称为内点记内点集为。边界与网格线的交点称为边界点,,h,边界点全体记为。与节点沿方向或方向只差一个步长的点(,)xy和(,)xy(,)xyy,xkj,kj,kjh,,,:称为节点(,)xy的相邻节点。如果一个内点的四个相邻节点均属于称为正则内点kj(),,:正内点的全体记为,至少有一个相邻节点不属于的内点称为非正则内点非正则内(),点的全体记为。问题是要求出第一边值问题在全体内点上的数值解。,,,,,,,,,,,,,,,,,,,,为简便记:(,)(,)kjxy,ukjuxy(,)(,),ffxy,(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,。对正则内点kjkjkjkj,()由二阶中心差商公式:(,)kj,,,,,uukjukjukjh(,)(,)(,)(),,,uxhy(,)kjx,xhkj(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uukjukjukj(,)(,)(,),(),,uxy(,),,kjy,y,kj(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(),,uuPoisson方程在点处可表示为:,fxy(,)(,)kj,,xyukjukjukjukjukjukj(,)(,)(,)(,)(,)(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,fRkj(,)kj,,h,,,,,,,,,,,()其中:h,()()(,)(,)(,),Rkjuxhyuxy,,,kjkjxx,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()(,)Oh,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(),,,,,,,,,,为其截断误差表示式略去即得与方程相近似的差分方程:R(k,j)uuuuuu,,kjkjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,fkj,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()u式中方程的个数等于正则内点的个数而未知数则除了包含正则内点处解的近似值外ukj,还包含一些非正则内点处的近似值因而方程个数少于未知数个数。在非正则内点处Poissonu方程的差分近似不能按上式给出需要利用边界条件得到。边界条件的处理可以有各种方案下面介绍较简单的两种。uu(),,,,,直接转移。用最接近非正则内点的边界点上的值作为该点上值的近似这就是边界条件的直接转移。例如点Pkj(,)为非正则内点其最接近的边界点为Q点则有:()uuQQkj,,,,()()(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,kj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()u上式可以看作是用零次插值得到非正则内点处的近似值容易求出其截断误差为。将上式代入方程个数即与未知数个数相等。O(h,)(),,,,,线性插值。这种方案是通过用同一条网格线上与点相邻的边界点与内点作线性插P值得到非正则内点处值的近似。由点R与T的线性插值确定的近似值得:uuPkj(,)uP()kj,hd,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uRTkj,hdhd,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()其中dRP,其截断误差为。将其与方程相近似的差分方程联立得到方程个数与未Oh()知数个数相等的方程组求解此方程组可得Poisson方程第一边值问题的数值解。上面所给出的差分格式称为五点菱形格式uuuuuu,,kjkjkjkjkjkj,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,fkj,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()h,,实际计算时经常取此时五点菱形格式可化为:uuuuuf,,(),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,kjkjkjkjkjkj,,,,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()简记为:uf,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,kjkj,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()uuuuuu,,其中:。kjkjkjkjkjkj,,,,,,,,例,,,,,,,,,,用五点菱形格式求解拉普拉斯(Laplace)方程第一边值问题。,uu,,fxyxy(,)(,),,,,xy,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uxyxy(,)lg(),,,,,,,其中:。取。,,,,,,,{(,),}xyxyh解:网格中有四个内点均为正则内点。由五点菱形格式得方程组:,()uuuuu,,,,,,,,h,,()uuuuu,,,,,,,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()uuuuu,,,,,,,,h,,()uuuuu,,,,,,,h,代入边界条件:,uu,,lg,lg,,,,,uu,,lg,lg,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uu,,lg,lg,,,,,uu,,lg,lg,,,其解为:u,u,u,u,,,,,h,,当时对利用点uuuuuf,,()(,)kj(,)kj,,(,)kj,kjkjkjkjkjkj,,,,,,,,h构造的差分格式称为五点矩形格式简记为:,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,uf,kjkj,,h,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()其中:其截断误差为:uuuuuu,,kjkjkjkjkjkj,,,,,,,,,,,,huuu,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,)()()RkjOhOh,,,,,,,,xxyy,,kj(,),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,()五点菱形格式与矩形格式的截断误差均为称它们具有二阶精度。如果用更多的点构Oh()造差分格式其截断误差的阶数可以提高如利用菱形格式及矩形格式所涉及的所有节点构造出的九点格式就是具有四阶精度的差分格式。,,,,,,,,,,用Dirichlet法求解Laplace方程的MATLAB程序,,,,,,,,,,,,,,,,,,,,求解区间uxtuxt(,)(,),内的近似解。而且满足Rxyxayb,,,,,{(,):,}xxyy,,xa条件:其中而且uxfxuxbfx(,)(),(,)(),,uyfyuayfy(,)(),(,)(),,其中而且存在整数和使得:。nm,,ybanhbmh,,,*****************************************************function,,,,,U=dirich(f,f,f,f,a,b,h,tol,max),,,,,Input,,,,,,,,,,f,f,f,f,,,,,are,,,,,boundary,,,,,functions,,,,,input,,,,,as,,,,,strings,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,and,,,,,b,,,,,right,,,,,endpoints,,,,,of,,,,,,a,,,,,and,,,,,,b,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,h,,,,,step,,,,,size,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,tol,,,,,is,,,,,the,,,,,toleranceOutput,,,,,,,,,,U,,,,,solution,,,,,matrix,,,,,analogous,,,,,to,,,,,Table,,,,,Initialize,,,,,parameters,,,,,and,,,,,U,,,,,n=fix(ah)m=fix(bh)ave=(a*(feval(f,)feval(f,)),,,,,,,,,,,,,,,,,,,,b*(feval(f,)feval(f,)))(*a*b)U=ave*ones(n,m),,,,,Boundary,,,,,conditions,,,,,U(,:m)=feval(f,:h:(m)*h)'U(n,:m)=feval(f,:h:(m)*h)'U(:n,)=feval(f,:h:(n)*h)U(:n,m)=feval(f,:h:(n)*h)U(,)=(U(,)U(,))U(,m)=(U(,m)U(,m))U(n,)=(U(n,)U(n,))U(n,m)=(U(n,m)U(n,m)),,,,,SOR,,,,,parameter,,,,,w=(sqrt((cos(pi(n))cos(pi(m)))^)),,,,,Refine,,,,,approximations,,,,,and,,,,,sweep,,,,,operator,,,,,throughout,,,,,the,,,,,grid,,,,,err=cnt=while((err>tol)(cnt<=max)),,,,,,,,,,,,,,,err=,,,,,,,,,,,,,,,for,,,,,j=:m,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,for,,,,,i=:n,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,relx=w*(U(i,j)U(i,j)U(i,j)U(i,j)*U(i,j)),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,U(i,j)=U(i,j)relx,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,if,,,,,(err<=abs(relx)),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,err=abs(relx),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,end,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,end,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,endcnt=cntend,,,,,U=flipud(U'

用户评价(0)

关闭

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

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

提示

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

文档小程序码

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

1

打开微信

2

扫描小程序码

3

发布寻找信息

4

等待寻找结果

我知道了
评分:

/45

Matlab解微分方程[整理版]

VIP

在线
客服

免费
邮箱

爱问共享资料服务号

扫描关注领取更多福利