机器学习计算方法常微分方程初值问题的数值解

计算方法第七章常微分方程初值问题地数值解7.1引入7.2常微分方程初值问题7.3欧拉方法及其改进7.4梯形方法7.5龙格-库塔方法7.6常微分方程组地数值解法2 27.1引入—三论PageRank算法设网页i地重要性为Pri,链出数为Li。如果网页i存在一个指向网页A地链接,则表明i地所有者认为A比较重要,从而把i地一部分重要性得分赋予A。一个页面地PageRank是由所有链向它地页面（链入页面）地重要性经过递归算法得到。假设世界上只有四张网页:A,B,C,D3 37.1引入—三论PageRank算法网页间地链接矩阵M为(7.1.1)其中M地元素mij为0时表示第i个网页没有到第j个网页地链接,否则mij为1/Li4 47.1引入—三论PageRank算法将网页地重要性Pr看作是不断变化地变量,则每次地变化为M·Pr-Pr=(M-I)·Pr。用建立下面地常微分方程地方法进行求解。(7.1.2)5 57.1引入—三论PageRank算法分别取Pr地初值为[0.25,0.25,0.25,0.25]T [1,0,0,0]T6 67.1引入常微分方程少数简单类型地常微分方程能求得精确解析有些常微分方程求解过程极为复杂多数情况只能使用近似解法求得其近似解常微分方程数值解法很有必要7 7第七章常微分方程初值问题地数值解7.1引入7.2常微分方程初值问题7.3欧拉方法及其改进7.4梯形方法7.5龙格-库塔方法7.6常微分方程组地数值解法8 87.2常微分方程初值问题一阶常微分方程地定解问题(7.2.1)其中:x∈[a,b]为自变量=(y1,y2,…,yd)∈Rd,y=y(x)为向量函数,(x,y):RRd→Rd称为右端向量场y0∈Rd称为初值。当给定向量场f(x,y)与初值y0,在区间[a,b]上求函数值y(x),使其满足方程组(7.2.1)地问题称为初值问题。在d=1地情况,方程组(7.2.1)为简单地常微分方程初值问题。9 9常微分方程初值问题地解地存在性与唯一性定理定理7.1假设f(x,y)在区域G={(x,y)|a≤x≤b,|y|<∞内连续,并对y满足,利普希茨(Lipschitz)条件,即存在常数L>0,使得(7.2.2)对所有x∈[a,b]与任意y1,y2∈R成立,则初值问题(7.2.1)在[a,b]上存在唯一解y(x),若f∈Ck(G),k≥0,那么其解y(x)∈Ck+1([a,b])107.2常微分方程初值问题例7.1马尔萨斯模型:马尔萨斯在分析人口出生与死亡情况地资料后发现,人口净增长率r基本上是一个常数,提出了著名地人口指数增长模型。(7.2.3)其中,t表示时间,N表示人口数117.2常微分方程初值问题例7.2一重物垂直作用于弹簧所引起地振荡,当运动阻力与速度平方成正比时,可用一个二阶常微分方程进行描述:(7.2.4)12例7.2若令则上述二阶常微分方程可以化为等价地一阶常微分方程组(7.2.5)13第七章常微分方程初值问题地数值解7.1引入7.2常微分方程初值问题7.3欧拉方法及其改进7.4梯形方法7.5龙格-库塔方法7.6常微分方程组地数值解法14 147.3欧拉方法及其改进对于初值问题(7.2.1),在求解区间[a,b]上取等间距节点a=x0<x1<x2<…<xN-1<xN=b(7.3.1)其中∆xn=xn+1-xn=h称为积分网格地步长常微分方程(7.2.1)初值问题地数值解采用数值算法计算出初值问题精确解y(x)在节点x1,x2,…,xN上地函数值地近似值y1,y2,…,yN。常用地数值算法:几何方法,数值微分,数值积分与Taylor展开15 157.3欧拉方法及其改进如果在计算yn时该迭代公式只用到已经求出地yn-1,而不使用y1,y2,…,yn-2中地任何一个则称此算法为单步方法否则称之为多步方法167.3.1欧拉方法17 177.3.1欧拉方法取等距节点组,此时有xn=a+nh,n=0,1,…,N,∆xn≡h=(b-a)/N若以向前差商近似代替(7.2.1)中地导数y‘(xn),得(7.3.2)187.3.1欧拉方法将其看作等式,并记y(xn)地近似值为yn,则有(7.3.3)19欧拉方法Octave程序functiony=Euler(f,y0,a,b,N)%EulermethodtosolveOrdinarydifferentialequationy'=f(x,y),a<x<=b%withinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,1);y(1)=y0+h*f(a,y0);fori=2:Ny(i)=y(i-1)+h*f(a+(i-1)*h,y(i-1));endendfunction207.3.2欧拉方法地几何解释21 217.3.2欧拉方法地几何解释前向欧拉方法地求解过程可以看作:从初值点(x0,y0)出发,以斜率f(x0,y0)做切线在经过步长h后与另一条积分曲线相交于(x1,y1)依次得到y1,y2,…,yN依次连接(x0,y0),(x1,y1),…,(xN,yN)地折线去近似方程y’(x)=f(x,y(x))过点(x0,y0)地积分曲线该方法又称为折线法227.3.2欧拉方法地几何解释237.3.3欧拉方法地截断误差24 247.3.3欧拉方法地截断误差假设y(x)是y’(x)=f(x,y(x))任意一个二次连续可微地解,则将y(x)在点xn处进行Taylor展开,可得(7.3.4)若取x=xn+1,则有(7.3.5)257.3.3欧拉方法地截断误差余项(7.3.6)为向前欧拉公式(7.3.3)地误差称为向前欧拉公式地截断误差267.3.4向后欧拉方法27 277.3.4向后欧拉方法用去近似导数有(7.3.7)记y(xn)地近似值为yn,则有(7.3.8)28"预报-校正"法先以向前欧拉公式给出yn+1地初值(7.3.9)然后再用迭代公式(7.3.10)计算yn+1地近似值。只用(7.3.10)式计算一次,而不进行迭代29"预报-校正"法Octave代码functiony=Euler(f,y0,a,b,N)EulermethodtosolveOrdinarydifferentialequationy'=f(x,y),a<x<=bwithinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,1);y(1)=y0+h*f(a,y0);fori=2:Ny(i)=y(i-1)+h*f(a+(i-1)*h,y(i-1));endendfunction30第七章常微分方程初值问题地数值解7.1引入7.2常微分方程初值问题7.3欧拉方法及其改进7.4梯形方法7.5龙格-库塔方法7.6常微分方程组地数值解法31 317.4.1梯形方法常微分方程地定解问题(7.2.1)式中地将其在[xn,xn+1]区间上进行积分,可得(7.4.1)如果等式右端地积分可以求出,就可以得到y(xn+1)可以采用不同地数值积分公式进行近似求解327.4.1梯形方法用梯形求积公式,可得(7.4.2)其中(7.4.3)为梯形求积公式余项337.4.1梯形方法记y(xn)地近似值为yn,则有(7.4.4)称此递推公式为梯形公式347.4.2改进欧拉格式35 357.4.2改进欧拉格式例7.3分别用向前欧拉法,基于"预报-校正"技术地向后欧拉法与改进欧拉格式求解初值问题36 36例7.3解:取步长h=0.1,则xn=0.1n(n=0,1,…,10),并由y0=1,分别采用三种方法进行求解。向前欧拉法地迭代公式为:37 37例7.3解:基于"预报-校正"技术地向后欧拉法地迭代公式为:38 38例7.3解:改进欧拉格式地迭代公式为:39 39例7.3地比较结果xn向前Euler法向后Euler法改进Euler解析解11.10001.09181.09591.095421.19181.17631.18411.183231.27741.25461.26621.264941.35821.32781.34341.341651.43511.39641.41641.414261.50901.46091.48601.483271.58031.52161.55251.549281.64981.57861.61651.612591.71781.63211.67821.6733101.78481.68191.73791.732140 40例7.3地比较结果41 41第七章常微分方程初值问题地数值解7.1引入7.2常微分方程初值问题7.3欧拉方法及其改进7.4梯形方法7.5龙格-库塔方法7.6常微分方程组地数值解法42 427.5.1龙格-库塔方法地基本思想43 437.5.1龙格-库塔方法地基本思想设y(x)是初值问题(7.2.1)地解,将y(xn+1)在xn处进行泰勒展开,只保留常数项,有(7.5.1)其中f(ξ,y(ξ))为y(x)在区间[xn,xn+1]上地平均斜率用区间[xn,xn+1]内更多节点地斜率地线性组合得到平均斜率地更高逼近44 447.5.2龙格-库塔方法45 457.5.2龙格-库塔方法设y(x)是初值问题(7.2.1)地解,并充分光滑,则将y(x)在xn处进行泰勒展开,取x=xn+1,有(7.5.2)其中(7.5.3)为余项467.5.2龙格-库塔方法将该余项略去,并记y(xn)地近似值为yn,则有(7.5.4)由微分地链式法则可知(7.5.5)477.5.2龙格-库塔方法将该余项略去,并记y(xn)地近似值为yn,则有(7.5.4)由微分地链式法则可知(7.5.5)487.5.2龙格-库塔方法用f(x,y)地函数值地线性组合近似(7.5.4)式右侧y(x)地各阶导数之与,即(7.5.6)其中(7.5.7)497.5.3二级龙格-库塔格式50 507.5.3二级龙格-库塔格式利用(7.5.5)式将(7.5.2)式改写为(7.5.8)其中f,f’x,f’y分别代表f(x,y(x))及其对x与y地一阶偏导数在(xn,yn)上取值。517.5.3二级龙格-库塔格式下面分析近似公式(7.5.6)与(7.5.7),级数M与参数αm,λm,μm地取值都不是唯一地不同地取法可以得到不同地迭代格式设M=2,由二元函数地泰勒展开式可以得到(7.5.9)527.5.3二级龙格-库塔格式将其代入(7.5.7)地第二式,进而将(7.5.7)式代入(7.5.6)式,可得(7.5.10)比较(7.5.8)与(7.5.10)两式中h与h2地系数,并令对应地各项系数相等,有(7.5.11)537.5.3二级龙格-库塔格式令则得二级龙格-库塔格式为(7.5.12)547.5.3二级龙格-库塔格式令则得二级Runge-Kutta格式为(7.5.13)55二级龙格-库塔格式Octave代码functiony=midpoint_RuggeKutta(f,y0,a,b,N)midpoint_RuggeKuttamethodtosolveODEy'=f(x,y),a<x<=bwithinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,1);y(1)=y0+h*f(a+h/2,y0+h*f(a,y0)/2);fori=2:Ny(i)=y(i-1)+h*(f(a+(i-1)*h+h/2,y(i-1))+h*f(a+(i-1)*h,y(i-1))/2);endendfunction567.5.4四级龙格-库塔格式57 577.5.4四级龙格-库塔格式当p=4时,取M=4(7.5.16)587.5.4四级龙格-库塔格式进一步利用Taylor展开式与(7.5.5)式可以推出(7.5.17)其中,是把(7.5.16)式ki中地yn与ki-1用y(xn)与代替后得到地表达式。并有(7.5.18)四级龙格-库塔格式地截断误差为h地5阶无穷小量59四级龙格-库塔格式Octave程序functiony=Fourth_RuggeKutta(f,y0,a,b,N)midpoint_RuggeKuttamethodtosolveODEy'=f(x,y),a<x<=bwithinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N+1,1);y(1)=y0;fori=2:N+1k1=f(a+(i-2)*h,y(i-1));k2=f(a+(i-2)*h+h/2,y(i-1)+k1*h/2);k3=f(a+(i-2)*h+h/2,y(i-1)+k2*h/2);k4=f(a+(i-1)*h,y(i-1)+k3*h);y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6;endendfunction607.5.4四级龙格-库塔格式例7.4用四阶龙格-库塔格式求解下面初值问题并与例7.3结果进行比较61例7.4解:取步长h=0.1,四阶龙格-库塔格式为62例7.4地结果比较xn向前Euler法向后Euler法改进Euler四阶Runge-解析解Kutta11.10001.09181.09591.09541.095421.19181.17631.18411.18321.183231.27741.25461.26621.26491.264941.35821.32781.34341.34161.341651.43511.39641.41641.41421.414261.50901.46091.48601.48321.483271.58031.52161.55251.54921.549281.64981.57861.61651.61251.612591.71781.63211.67821.67331.6733101.78481.68191.73791.73211.732163第七章常微分方程初值问题地数值解7.1引入7.2常微分方程初值问题7.3欧拉方法及其改进7.4梯形方法7.5龙格-库塔方法7.6常微分方程组地数值解法64 647.6常微分方程组地数值解法65 657.6常微分方程组地数值解法一阶常微分方程组(7.6.1)其向前欧拉格式为(7.6.2)667.6常微分方程组地数值解法改进地欧拉格式为(7.6.3)四阶龙格-库塔格式为(7.6.4)677.6常微分方程组地数值解法例7.5采用向前欧拉格式,