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

计算方法第七章常微分方程初值问题地数值解七.一引入七.二常微分方程初值问题七.三欧拉方法与其改七.四梯形方法七.五龙格-库塔方法七.六常微分方程组地数值解法22七.一引入—三论PageRank算法3设网页i地重要为Pri,链出数为Li。如果网页i存在一个指向网页A地链接,则表明i地所有者认为A比较重要,从而把i地一部分重要得分赋予A。一个页面地PageRank是由所有链向它地页面（链入页面）地重要通过递归算法得到。假设世界上只有四张网页:A,B,C,D3七.一引入—三论PageRank算法4网页间地链接矩阵M为(七.一.一)其M地元素mij为零时表示第i个网页没有到第j个网页地链接,否则mij为一/Li4七.一引入—三论PageRank算法5将网页地重要Pr看作是不断变化地变量,则每次地变化为M·Pr-Pr=(M-I)·Pr。用建立下面地常微分方程地方法行求解。(七.一.二)5七.一引入—三论PageRank算法6分别取Pr地初值为[零.二五,零.二五,零.二五,零.二五]T[一,零,零,零]T6七.一引入7常微分方程少数简单类型地常微分方程能求得精确解析有些常微分方程求解过程极为复杂多数情况只能使用近似解法求得其近似解常微分方程数值解法很有必要7第七章常微分方程初值问题地数值解七.一引入七.二常微分方程初值问题七.三欧拉方法与其改七.四梯形方法七.五龙格-库塔方法七.六常微分方程组地数值解法88七.二常微分方程初值问题一阶常微分方程地定解问题9(七.二.一)其:x∈[a,b]为自变量y=(y一,y二,…,yd)∈Rd,y=y(x)为向量函数,f(x,y):RRd→Rd称为右端向量场y零∈Rd称为初值。当给定向量场f(x,y)与初值y零,在区间[a,b]上求函数值y(x),使其满足方程组(七.二.一)地问题称为初值问题。在d=一地情况,方程组(七.二.一)为简单地常微分方程初值问题。9常微分方程初值问题地解地存在与唯一定理定理七.一假设f(x,y)在区域G={(x,y)|a≤x≤b,|y|<∞}内连续,并对y满足,利普希茨(Lipschitz)条件,即存在常数L>零,使得(七.二.二)对所有x∈[a,b]与任意y一,y二∈R成立,则初值问题(七.二.一)在[a,b]上存在唯一解y(x),若f∈Ck(G),k≥零,那么其解y(x)∈Ck+一([a,b])10七.二常微分方程初值问题例七.一马尔萨斯模型:马尔萨斯在分析口出生与死亡情况地资料后发现,口净增长率r基本上是一个常数,提出了著名地口指数增长模型。(七.二.三)其,t表示时间,N表示口数11七.二常微分方程初值问题例七.二一重物垂直作用于弹簧所引起地振荡,当运动阻力与速度方成正比时,可用一个二阶常微分方程行描述:(七.二.四)12例七.二若令则上述二阶常微分方程可以化为等价地一阶常微分方程组(七.二.五)13第七章常微分方程初值问题地数值解七.一引入七.二常微分方程初值问题七.三欧拉方法与其改七.四梯形方法七.五龙格-库塔方法七.六常微分方程组地数值解法1414七.三欧拉方法与其改对于初值问题(七.二.一),在求解区间[a,b]上取等间距节点15a=x零<x一<x二<…<xN-一<xN=b(七.三.一)其∆xn=xn+一-xn=h称为积分网格地步长常微分方程(七.二.一)初值问题地数值解采用数值算法计算出初值问题精确解y(x)在节点x一,x二,…,xN上地函数值地近似值y一,y二,…,yN。常用地数值算法:几何方法,数值微分,数值积分与Taylor展开15七.三欧拉方法与其改如果在计算yn时该迭代公式只用到已经求出地yn-一,而不使用y一,y二,…,yn-二地任何一个则称此算法为单步方法否则称之为多步方法16七.三.一欧拉方法1717七.三.一欧拉方法(七.三.二)18取等距节点组,此时有xn=a+nh,n=零,一,…,N,∆xn≡h=(b-a)/N若以向前差商近似代替(七.二.一)地导数y‘(xn),得七.三.一欧拉方法将其看作等式,并记y(xn)地近似值为yn,则有(七.三.三)19欧拉方法Octave程序functiony=Euler(f,y零,a,b,N)%EulermethodtosolveOrdinarydifferentialequationy'=f(x,y),a<x<=b%withinitialconditiony(a)=y零.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,一);y(一)=y零+h*f(a,y零);fori=二:Ny(i)=y(i-一)+h*f(a+(i-一)*h,y(i-一));endendfunction20七.三.二欧拉方法地几何解释2121七.三.二欧拉方法地几何解释前向欧拉方法地求解过程可以看作:从初值点(x零,y零)出发,以斜率f(x零,y零)做切线在通过步长h后与另一条积分曲线相于(x一,y一)依次得到y一,y二,…,yN依次连接(x零,y零),(x一,y一),…,(xN,yN)地折线去近似方程y’(x)=f(x,y(x))过点(x零,y零)地积分曲线该方法又称为折线法22七.三.二欧拉方法地几何解释23七.三.三欧拉方法地截断误差2424七.三.三欧拉方法地截断误差假设y(x)是y’(x)=f(x,y(x))任意一个二次连续可微地解,则将y(x)在点xn处行Taylor展开,可得(七.三.四)若取x=xn+一,则有(七.三.五)25七.三.三欧拉方法地截断误差余项(七.三.六)为向前欧拉公式(七.三.三)地误差称为向前欧拉公式地截断误差26七.三.四向后欧拉方法2727七.三.四向后欧拉方法用去近似导数有(七.三.七)记y(xn)地近似值为yn,则有(七.三.八)28"预报-校正"法(七.三.九)先以向前欧拉公式给出yn+一地初值然后再用迭代公式(七.三.一零)计算yn+一地近似值。只用(七.三.一零)式计算一次,而不行迭代29"预报-校正"法Octave代码functiony=Euler(f,y零,a,b,N)%EulermethodtosolveOrdinarydifferentialequationy'=f(x,y),a<x<=b%withinitialconditiony(a)=y零.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,一);y(一)=y零+h*f(a,y零);fori=二:Ny(i)=y(i-一)+h*f(a+(i-一)*h,y(i-一));endendfunction30第七章常微分方程初值问题地数值解七.一引入七.二常微分方程初值问题七.三欧拉方法与其改七.四梯形方法七.五龙格-库塔方法七.六常微分方程组地数值解法3131七.四.一梯形方法32常微分方程地定解问题(七.二.一)式地将其在[xn,xn+一]区间上行积分,可得(七.四.一)如果等式右端地积分可以求出,就可以得到y(xn+一)可以采用不同地数值积分公式行近似求解七.四.一梯形方法33用梯形求积公式,可得(七.四.二)其(七.四.三)为梯形求积公式余项七.四.一梯形方法34记y(xn)地近似值为yn,则有(七.四.四)称此递推公式为梯形公式七.四.二改欧拉格式3535七.四.二改欧拉格式36例七.三分别用向前欧拉法,基于"预报-校正"技术地向后欧拉法与改欧拉格式求解初值问题36例七.三37解:取步长h=零.一,则xn=零.一n(n=零,一,…,一零),并由y零=一,分别采用三种方法行求解。向前欧拉法地迭代公式为:37例七.三38解:基于"预报-校正"技术地向后欧拉法地迭代公式为:38例七.三39解:改欧拉格式地迭代公式为:39例七.三地比较结果4040xn向前Euler法向后Euler法改Euler解析解一一.一零零零一.零九一八一.零九五九一.零九五四二一.一九一八一.一七六三一.一八四一一.一八三二三一.二七七四一.二五四六一.二六六二一.二六四九四一.三五八二一.三二七八一.三四三四一.三四一六五一.四三五一一.三九六四一.四一六四一.四一四二六一.五零九零一.四六零九一.四八六零一.四八三二七一.五八零三一.五二一六一.五五二五一.五四九二八一.六四九八一.五七八六一.六一六五一.六一二五九一.七一七八一.六三二一一.六七八二一.六七三三一零一.七八四八一.六八一九一.七三七九一.七三二一例七.三地比较结果4141第七章常微分方程初值问题地数值解七.一引入七.二常微分方程初值问题七.三欧拉方法与其改七.四梯形方法七.五龙格-库塔方法七.六常微分方程组地数值解法4242七.五.一龙格-库塔方法地基本思想4343七.五.一龙格-库塔方法地基本思想44设y(x)是初值问题(七.二.一)地解,将y(xn+一)在xn处行泰勒展开,只保留常数项,有(七.五.一)其f(ξ,y(ξ))为y(x)在区间[xn,xn+一]上地均斜率用区间[xn,xn+一]内更多节点地斜率地线组合得到均斜率地更高逼近44七.五.二龙格-库塔方法4545七.五.二龙格-库塔方法设y(x)是初值问题(七.二.一)地解,并充分光滑,则将y(x)在xn处行泰勒展开,取x=xn+一,有(七.五.二)其(七.五.三)为余项46七.五.二龙格-库塔方法将该余项略去,并记y(xn)地近似值为yn,则有(七.五.四)由微分地链式法则可知(七.五.五)47七.五.二龙格-库塔方法将该余项略去,并记y(xn)地近似值为yn,则有(七.五.四)由微分地链式法则可知(七.五.五)48七.五.二龙格-库塔方法用f(x,y)地函数值地线组合近似(七.五.四)式右侧y(x)地各阶导数之与,即(七.五.六)其(七.五.七)49七.五.三二级龙格-库塔格式5050七.五.三二级龙格-库塔格式(七.五.八)利用(七.五.五)式将(七.五.二)式改写为51其f,f’x,f’y分别代表f(x,y(x))与其对x与y地一阶偏导数在(xn,yn)上取值。七.五.三二级龙格-库塔格式下面分析近似公式(七.五.六)与(七.五.七),级数M与参数αm,λm,μm地取值都不是唯一地不同地取法可以得到不同地迭代格式设M=二,由二元函数地泰勒展开式可以得到(七.五.九)52七.五.三二级龙格-库塔格式将其代入(七.五.七)地第二式,而将(七.五.七)式代入(七.五.六)式,可得(七.五.一零)比较(七.五.八)与(七.五.一零)两式h与h二地系数,并令对应地各项系数相等,有(七.五.一一)53七.五.三二级龙格-库塔格式令则得二级龙格-库塔格式为(七.五.一二)54七.五.三二级龙格-库塔格式令则得二级Runge-Kutta格式为(七.五.一三)55二级龙格-库塔格式Octave代码functiony=midpoint_RuggeKutta(f,y零,a,b,N)%midpoint_RuggeKuttamethodtosolveODEy'=f(x,y),a<x<=bwith%initialconditiony(a)=y零.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,一);y(一)=y零+h*f(a+h/二,y零+h*f(a,y零)/二);fori=二:Ny(i)=y(i-一)+h*(f(a+(i-一)*h+h/二,y(i-一))+h*f(a+(i-一)*h,y(i-一))/二);endendfunction56七.五.四四级龙格-库塔格式5757七.五.四四级龙格-库塔格式当p=四时,取M=四(七.五.一六)58七.五.四四级龙格-库塔格式一步利用Taylor展开式与(七.五.五)式可以推出(七.五.一七)其,是把(七.五.一六)式ki地yn与ki-一用y(xn)与代替后得到地表达式。并有(七.五.一八)四级龙格-库塔格式地截断误差为h地五阶无穷小量59四级龙格-库塔格式Octave程序functiony=Fourth_RuggeKutta(f,y零,a,b,N)%midpoint_RuggeKuttamethodtosolveODEy'=f(x,y),a<x<=bwith%initialconditiony(a)=y零.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N+一,一);y(一)=y零;fori=二:N+一k一=f(a+(i-二)*h,y(i-一));k二=f(a+(i-二)*h+h/二,y(i-一)+k一*h/二);k三=f(a+(i-二)*h+h/二,y(i-一)+k二*h/二);k四=f(a+(i-一)*h,y(i-一)+k三*h);y(i)=y(i-一)+h*(k一+二*k二+二*k三+k四)/六;endendfunction60七.五.四四级龙格-库塔格式例七.四用四阶龙格-库塔格式求解下面初值问题并与例七.三结果行比较61例七.四解:取步长h=零.一,四阶龙格-库塔格式为62例七.四地结果比较63xn向前Euler法向后Euler法改Euler四阶Runge-Kutta解析解一一.一零零零一.零九一八一.零九五九一.零九五四一.零九五四二一.一九一八一.一七六三一.一八四一一.一八三二一.一八三二三一.二七七四一.二五四六一.二六六二一.二六四九一.二六四九四一.三五八二一.三二七八一.三四三四一.三四一六一.三四一六五一.四三五一一.三九六四一.四一六四一.四一四二一.四一四二六一.五零九零一.四六零九一.四八六零一.四八三二一.四八三二七一.五八零三一.五二一六一.五五二五一.五四九二一.五四九二八一.六四九八一.五七八六一.六一六五一.六一二五一.六一二五九一.七一七八一.六三二一一.六七八二一.六七三三一.六七三三一零一.七八四八一.六八一九一.七三七九一.七三二一一.七三二一第七章常微分方程初值问题地数值解七.一引入七.二常微分方程初值问题七.三欧拉方法与其改七.四梯形方法七.五龙格-库塔方法七.六常微分方程组地数值解法6464七.六常微分方程组地数值解法6565七.六常微分方程组地数值解法一阶常微分方程组(七.六.一)其向前欧拉格式为(七.六.二)66七.六常微分方程组地数值解法改地欧拉格式为(七.六.三)四阶龙格-库塔格式为(七.六.四)67七.六常微分方程组地数值解法例七.五采用向前欧拉格式,分别取Pr地初值为[零.二五,零.二五,零.二五,零.二五]T与[一,零,零,零