流体传热与数值模拟第九讲

流体流动与传热过程数值模拟讲义之九/*NumericalMethodsforOrdinaryDifferentialEquations*/常微分方程数值解主要内容•常微分方程概念•数值求解方法–隐式法–显式法•欧拉法（重点）•龙格-库塔法（重点2）•应用举例–欧拉及R-K法编程求解例：对下图所示的二维非稳态金属块的无源导热问题，其初始温度为室温27度，四周边界温度为定值。求解金属块内的非稳态二维热传导。导热系数k=1000w/(mK)。控制方程：dydTdxddxdTdxdcktT边界条件：0x100T5.0x500Tx0y200T5.0y300TY300NT200ST100WT500ET常微分方程提出初始条件：t=0,T=2721,,1,2,1,,1,,,22)(yTTTxTTTadydTkdxddxdTkdxdaTfjijijijijijijijijidydTkdxddxdTkdxdaTf)(dydTkdxddxdTkdxdatT)(TftT定义得到常微分方程（对时间微分）（对空间微分）初值条件t=0,T=27常微分方程与空间微分方程区别•对时间或演变过程因变量进行微分-推进或演变概念。•初值与边界条件常微分方程数值解/*NumericalMethodsforOrdinaryDifferentialEquations*/待求解的问题：一阶常微分方程的初值问题/*Initial-ValueProblem*/:0)(],[),(yaybaxyxfdxdy解的存在唯一性（“常微分方程”理论）：只要f(x,y)在[a,b]R1上连续，且关于y满足Lipschitz条件，即存在与x,y无关的常数L使对任意定义在[a,b]上的y1(x)和y2(x)都成立，则上述IVP存在唯一解。|||),(),(|2121yyLyxfyxf解析解法：（常微分方程理论）只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。如何求解计算解函数y(x)在一系列节点a=x0x1…xn=b处的近似值),...,1()(nixyyii节点间距为步长，通常采用等距节点，即取hi=h(常数)。)1,...,0(1nixxhiii数值解法:求解所有的常微分方程步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。111()()()()()(,)nnnnnnnnnnyxyxhyxyxyyxyyhfxy1(,)0,1,...nnnnyyhfxyn--------Euler’sMethod§1欧拉方法/*Euler’sMethod*/§1Euler’sMethodTaylor展开法几何意义亦称为欧拉折线法/*Euler’spolygonalarcmethod*/几何直观是帮助我们寻找解决一个问题的思路的好办法哦定义在假设yn=y(xn)，即第n步计算是精确的前提下，考虑公式或方法本身带来的误差:Rn=y(xn+1)yn+1,称为局部截断误差/*localtruncationerror*/。说明显然，这种近似有一定误差，而且步长越大，误差越大，如何估计这种误差y(xn+1)yn+1？§1Euler’sMethod截断误差:实际上，y(xn)yn，yn也有误差，它对yn+1的误差也有影响，见下图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。局部截断误差的分析：由于假设yn=y(xn)，即yn准确，因此分析局部截断误差时将y(xn+1)和yn+1都用点xn上的信息来表示，工具：Taylor展开。欧拉法的局部截断误差：223111232()[()()()()][(,)]()()hnnnnnnnnnhnRyxyyxhyxyxOhyhfxyyxOhRn+1的主项/*leadingterm*/§1Euler’sMethod定义若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。欧拉法具有1阶精度。在xn点用一阶向前差商近似一阶导数1()()()nnnyxyxyxh在讨论差分公式时介绍了差商的概念和性质，各阶差商可以近似各阶导数，具有不同的精度，且可以用函数值来表示。111()()()()()(,)nnnnnnnnnnyxyxhyxyxyyxyyhfxyEuler’smethod§1Euler’sMethod§1Euler’sMethod欧拉公式的改进：隐式欧拉法或后退欧拉法/*implicitEulermethodorbackwardEulermethod*/11()()()nnnyxyxyxhxn+1点向后差商近似导数111111()()()()()(,)nnnnnnnnnnyxyxhyxyxyyxyyhfxy隐式或后退欧拉公式由于未知数yn+1同时出现在等式的两边，故称为隐式/*implicit*/欧拉公式，而前者称为显式/*explicit*/欧拉公式。隐式公式不能直接求解，一般需要用Euler显式公式得到初值，然后用Euler隐式公式迭代求解。因此隐式公式较显式公式计算复杂，但稳定性好（后面分析）。01(1)()111(,)(,)nnnnkknnnnyyhfxyyyhfxy(1)()1111111()(0)1111(1)11111()1(,)(,)1,()(,)kknnnnnnkknnnnknnnnnnknyyhfxyfxyhLyyhLyyhLyykyyhfxyy在迭代公式中取极限，有因此的极限就是隐式方程的解收敛性§1Euler’sMethod见上图，显然，这种近似也有一定误差，如何估计这种误差y(xn+1)yn+1？方法同上，基于Taylor展开估计局部截断误差。但是注意，隐式公式中右边含有f(xn+1,yn+1),由于yn+1不准确，所以不能直接用y'(xn+1)代替f(xn+1,yn+1)设已知曲线上一点Pn(xn,yn),过该点作弦线，斜率为(xn+1,yn+1)点的方向场f(x,y)方向,若步长h充分小，可用弦线和垂线x=xn+1的交点近似曲线与垂线的交点。几何意义xnxn+1PnPn+1xyy(x)§1Euler’sMethod1(,)0,1,...nnnnyyhfxyn比较尤拉显式公式和隐式公式及其局部截断误差231112()()()hnnnnRyxyyxOh显式公式111(,)nnnnyyhfxy隐式公式231112()()()hnnnnRyxyyxOh§1Euler’sMethod若将这两种方法进行算术平均，即可消除误差的主要部分/*leadingterm*/而获得更高的精度,称为梯形法§1Euler’sMethod梯形公式/*trapezoidformula*/—显、隐式两种算法的平均111[(,)(,)]2nnnnnnhyyfxyfxy注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。3111()()nnnRyxyOh梯形法的迭代计算和收敛性01(1)()111(,)(,)(,)2nnnnkknnnnnnyyhfxyhyyfxyfxy(1)()1111111()(0)1111(1)11111()1(,)(,)22221,()2(,)(,)2kknnnnnnkknnnnknnnnnnnnknhyyfxyfxyhhLLyyyyhLhyykLhyyfxyfxyy当时，在迭代公式中取极限，有因此的极限就是隐式方程的解收敛性§1Euler’sMethod梯形法的简化计算迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下预测校正系统，也称为改进尤拉法：改进欧拉法/*modifiedEuler’smethod*/Step1:先用显式欧拉公式作预测，算出),(1nnnnyxfhyyStep2:再将代入隐式梯形公式的右边作校正，得到1ny)],(),([2111nnnnnnyxfyxfhyy11(,),(,)2nnnnnnnnhyyfxyfxyhfxy§1Euler’sMethod注：此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。§1Euler’sMethod+1+11+1(,)(,(,))21(,)(,))2nnnnnnnnnnnnnnnhyyfxyfxhyhfxyyyhfxyyhfxy或其它形式1+1(,)(,)12pnnncnnpnpcyyhfxyyyhfxyyyy或几何解释xnxn+1ABPn+1=(A+B)/2尤拉法后退尤拉法梯形法§1Euler’sMethod中点欧拉公式/*midpointformula*/中心差商近似导数hxyxyxy2)()()(021x0x2x1))(,(2)()(1102xyxfhxyxy1,...,1),(211niyxfhyyiiii假设，则可以导出即中点公式也具有2阶精度，且是显式的。)(),(11iiiixyyxyy)()(311hOyxyRiii需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。§1Euler’sMethod几何解释xnxn+1尤拉法后退尤拉法中点法xn-10)(],[),(yaybaxyxfdxdy00(),xxyxyftytdt令x=x2,得2020(),xxyxyftytdtAnotherpointofview对右端积分采用中矩形公式即得中点公式§1Euler’sMethod预测-校正-改进系统中点法具有二阶精度，且是显式的，与梯形公式精度相匹配，用中点公式作预测，梯形公式作校正，得到如下预测校正系统：111112(,)(,)(,)2nnnnnnnnnnyyhfxyhyyfxyfxy预测：校正：3(3)1113(3)111312nnnnnnnnnnhyyyxyyxhyyyxyyx预测误差(设,准确):校正误差(设,准确):111114nnnnyxyyxy校正误差约为预测误差的1/4§1Euler’sMethod11111111111454515nnnnnnnnnnnyxyyyxyyyyxyyy预测误差和校正误差的事后误差估计式利用上两式可以估计预测值和校正值与准确值的误差，可以期望，利用这两个误差分别作预测值和校正值的补偿，有可能提高精度。设pn,cn分别为第n步的预测