（李荣荣）129基于优化方法的常微分方程边值问题数值解

（李荣荣）129基于优化方法的常微分方程边值问题数值解 12.9 基于优化方法的常微分方程边值问题数值解 常微分方程边值问题的一般形式为: ,, (12-2) yfxy,(,) (12-3) BCyaybxab((),())0,,,,,, BC式中: 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示边界条件所满足的函数关系。 12.9.1 基于Matlab函数的求解方法 Matlab求解边值问题的函数为bvp4c，它采用有限差分法求解，其基本 格式 pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载 为: solinit = bvpinit(x, yinit, params) y = bvp4c(odefun,bcfun,solinit) 函数bvpinit输入参数依次为自变量的区间[],函数的一个猜测值。 yxab, 函数bvp4c的输入参数依次为一阶微分方程或一阶微分方程组，用函数odefun定义， 边界条件用函数bcfun定义，这两个函数名用户自行定义。 【例12-1】 求解二阶常微分方程 ,,yy，,0 y(0)0, y(4)2,, 解:令 yy,1 ,yy, 12 ,yy,,21 yy的取值区间为[]，和的猜测值分别为0和2。计算程序为: x0,412 function ode_bvp1 clc;clf;clear all; solinit = bvpinit(linspace(0,4,5),[1 0]); sol = bvp4c(@twoode,@twobc,solinit); x = linspace(0,4); y = deval(sol,x); plot(x,y(1,:));%图12.18 function dydx = twoode(x,y) dydx = [ y(2) -abs(y(1))]; 1 function res = twobc(ya,yb) res = [ ya(1) yb(1)+ 2]; 图12.18 例12.1 数值解曲线 12.9.2 求解两点边值问题的打靶法 打靶法采用求初值的方法求解边值问题。用求解初值的方法求解式(12-2)和式(12-3) ,,描述的问题时，缺少处的条件。如果能找到合适的值，通过求解恰好能xa,ya()ya() *,使计算出的等于给定的边界值，则问题就得到解决。显然，寻找的过程yb()ya()yb() 是一个反复迭代的过程。就好比打靶，子弹以合适的角度和初速度射出才能命中目标。这就是打靶法的基本思想。 ,,,设一两点边值问题为: ， (12-4) yfxyy,(,,)yayb()(),,,,， ,,,,对应的初值问题为: ， (12-5) yfxyy,(,,)yayau()(),,,， 由该初值问题求解得出的另一端边界值为:ybu()(),,; u则ruu()()0,,,,,应满足的方程为: 。 因此，应用求初值问题的方法求边值问题转换成求上述方程的根的问题。求一元函数 2 根的方法很多，下面是用插值修正的方法求初始值，其迭代格式为: u ,,,()u2 uuuu,，,()2221()()uu,,,21 ,,,【例12-2】 求二阶常微分方程的解。 yyyyy，,,,30,(0)0,(2)1 解:该方程对应的一阶微分方程组为: ,yy，，，，12,y,, ,,,,,yyy,3212，，，， ,初始条件为:(猜测值)，计算程序如下 yy(0)0,(0)1,, function shoot1 clc; ii=0; tol=1e-9; xspan=[0,2]; guess=1; g1=guess; target=1; y0=[0;g1]; [x,y]=ode45(@dEqs,xspan,y0); t1=y(end,1) subplot(2,1,1); plot(x,y(:,1)) g2=1.1*g1; y0=[0;g2]; [x,y]=ode45(@dEqs,xspan,y0); t2=y(end,1) subplot(2,1,2); plot(x,y(:,1)) while abs(t2-target)>tol g=g2+(g2-g1)*(target-t2)/(t2-t1); t1=t2; y0=[0,g]; [x,y]=ode45(@dEqs,xspan,y0); t2=y(end,1); g1=g2; 3 g2=g; g,g1,g2,t1,t2 ii=ii+1; if ii>200 break end end figure(2) plot(x,y(:,1)) x,y function F=dEqs(x,y) % First-order differential F=[y(2);-3*y(1)*y(2)]; % equations. 计算结果为: ,,y(0)=0,1.5145,y(2)=1, 0.0145。 y(0),y(2), 12.9.3 边界层微分方程组及相似解 边界层微分方程组包括连续性方程，动量方程和能量方程，在二维稳态流动情况下: 对于强迫对流 ,,uv连续性方程为 (12-6) ，,0 ,,xy 2,,,uudpu1动量方程为 (12-7) ,u，,,，v2,,,,xydxy 12puconst，,,在边界层外缘 (12-8) ss2 2tt,,,,,,,w111,能量方程为 ， 无量纲温度 (12-9) ,ua，,v12tt,xyy,,,sw 对于自然对流，动量方程为 2,,,uudpu1 (12-10) ,ug，,,，,v2,,,,xydxy dp,,,g在边界层外缘 (12-11) sdx 4 tt,s能量方程的形式与(12-9)式相同，而无量纲温度通常表示为:。 ,,2tt,ws通过引入流函数和相似变换，动量方程和能量方程变为只依赖单个相似变量的常微分 方程。 ,,,,对于强迫对流 (12-12) u,,,,v yx,, ,f(),无量纲流函数 (12-13) , ux,s y,,相似变量 (12-14) xu/,s m取外流速度，动量方程变为 ucx,s m，12,,,,,, (12-15) fffmf，，,,(1)0 2 m,0此式是Falkner-Skan方程的另一种形式，当时，则变成Blasius方程。式中， ,m,为楔形体夹角。 ,, 2,, 1,,,,,，，,能量方程变为 (1)Pr0mf (12-16) 112 ,,,,,,,0,(0)(0)(0)0ff,1边界条件为 (12-17) ,,,,,,,,,()()1f,,1, 对于自然对流 1 4Gr,,x()/([4]),,,无量纲流函数 f, (12-18) ,,4,, 1 4Gry,,x,,相似变量 (12-19) ,,4x,, 2,,,,,,ffff320,动量方程变为 ，,，, (12-20) 2 5 ,,,能量方程变为 (12-21) ,,3Pr0，,f22 ,,,,,,,0,(0)(0)0,(0)1ff,2边界条件为 (12-22) ,,,,,,,,,()()0f,,2, 12.9.4 流函数方程和温度方程的求解 1. 应用打靶法求解Blasius方程 方程(12-15)、(12-16)、(12-20)、(12-21)为带渐近边界条件的两点边值问题，可采 用级数法、差分法、微分方程降阶法等方法求解。用得较多，且精度较高的方法是将高阶 微分方程化成一阶微分方程组，再用合适的微分方程数值方法求解。这样处理后，原两点 边值问题变成带未知初始条件的初值问题，需用试探或打靶的方法求解。 以强迫对流为例，初始条件为 ,,, = 未知 ,0,(0)(0)0,(0)fff,,, 未知初始条件的选取，应满足相应的远端边界条件。 应用打靶法求解Blasius方程的Matlab程序 %solution of Blasius function by shooting method %blasius_shoot.m clc; clear all; close all; tol=1e-9; xspan=[0,10]; guess=1; g1=guess; target=1; y0=[0;0;g1]; [x,y]=ode45('blasius',xspan,y0,[]); t1=y(end,2) subplot(2,1,1); plot(x,y(:,2)) g2=1.1*g1; y0=[0;0;g2]; [x,y]=ode45('blasius',xspan,y0,[]); t2=y(end,2) subplot(2,1,2); plot(x,y(:,2)) while abs(t2-target)>tol 6 g=g2+(g2-g1)*(target-t2)/(t2-t1); t1=t2; y0=[0,0,g]; [x,y]=ode45('blasius',xspan,y0,[]); t2=y(end,2); g1=g2; g2=g; g,g1,g2,t1,t2 i=i+1; if i>20 break end end figure(2) plot(x,y(:,2)) y 结果: y=[ 0 0 0.3320] 2. 实现打靶法的优化算法 寻找满足方程(12-15)、(12-16)、(12-20)、(12-21)远端边界条件的初始值的 过程可以看作是优化设计问题，即 ,,,求 使 xf,(0),(0),,, ,,, (12-23) min((0),(0))Ff, 目标函数由远端边界条件构成，可以取如下几种形式。 对强迫对流 ,,,, (12-24) Fff((0),(0))()1()1,,,,,，,,1 22,,,, (12-25) Fff((0),(0))(()1)(()1),,,,,，,,1 2222,,,,,,, (12-26) Ffff((0),(0))(()1)(())(()1)(()),,,,,,，,，,,，,11对自然对流边界层微分方程，目标函数构成与式(12-23)~(12-26)类似。 3. 与边界层微分方程对应的一阶微分方程组及目标函数 function ff=blasius(x,y) ff=[y(2);y(3);-0.5*y(1)*y(3)]; function f=blasiust(x,y,pr) 7 f=[y(2);y(3);-0.5*y(1)*y(3);y(5);-0.5*pr*y(1)*y(5)]; function fn=blasiust2(x,pr,eta_max) xspan=[0 eta_max]; y0=[0 0 x(1) 0 x(2)]; [eta ff]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); fn=(1-ff(end,2))^2+(1-ff(end,4))^2+(ff(end,3)^2+(ff(end,5)))^2; function f=naturalt(x,y,flag,pr) f=[y(2);y(3);-3*y(1)*y(3)+2*y(2)^2-y(4);y(5);-3*pr*y(1)*y(5)]; function fn=naturalt21(x,pr,eta_max) xspan=[0 eta_max]; y0=[0 0 x(1) 1 x(2)]; [eta ff]=ode45('naturalt',xspan,y0,[],pr); fn=[(ff(end,2))^2+(ff(end,3))^2+(ff(end,4))^2+(ff(end,5))^2]; 4. 应用纳齐谢姆-斯威格特方法求解 1)纳齐谢姆-斯威格特(Nachtsheim-Swigert)方法 纳齐谢姆-斯威格特法通过最小二乘法求出修正初值的增量，它既是一种实现打靶的方 法，也可看成是一种优化算法。 ,,,,,,记，，，。这样远端边界值可看成初y,(0),,,xf,(0)y,,(0),,,xf(0)1 ,,,始值的函数，即:，，，ffxy()(,),,ffxy()(,),,,()(,),,fxy1213 ,。 ,()(,),,fxy14 根据式(12-26)，由假设的初始值产生的远端边界值的误差可表示为: ,,(,)()1xyf,,,,1,,,(,)()xyf,,,,2 (12-27) ,(,)()1xy,,,,,31, ,,(,)()xy,,41,,, 2222误差的平方和为: (12-28) E,，，，,,,,1234 使式(12-28)为最小的初值就是满足远端边界条件的初值。根据函数取极值的必要条 件可得: ,,,,,,,,,31240，，，,,,,,1234,xxxx,,,,, (12-29) ,,,,,,,,,3124,0，，，,,,,,1234yyyy,,,,,, 8 通过对式(12-27)中的远端边界条件作一阶泰勒展开，并带入式(12-29)可得初值修正增量为: QQQQ,,,,yycxxycy,,x,2QQQ,,xxyyxy, (12-30) ,QQQQ,,,xxcyxycx,,,y2,QQQ,,xxyyxy, 式中 2222,,,,,，，，，Qff,,，,，,，,()()()(),,,,,,xxxxxx1,1,，，，，, ,,,,,,,Qffff,,,,，,,,，,,,()()()()()(),,xyxyxyxy1,1,,,,，,,,()(),,,1,1,xy ,,,,,,,Qffff,,,,,,,,,，,,,(1())()()()(1()),,cxxx1 (12-31) ,,,()()(),,,,,,1,11,xx,,,, 2222,,,,,，，，，，，，，f()()(),,Qf,,，(),，,，,yyy1,1,yyy,，，，，，，，， ,,,,,,,Qffff(1())()()()(1()),,,,,,,,,，,,,,cyyy1, ,,,()()(),,,,,,,,1,11,yy, 初始值的迭代格式为: (1)()()kkk，,,,,,,,fff(0)(0)(0),，,,,,,,,, (12-32) ,(1)()()kkk，,,,,,,(0)(0)(0),，,,,,,,,,111, 2)纳齐谢姆-斯威格特方法的Matlab程序 用打靶法求解式(12-15)、式(12-16)时需反复调用解常微分方程初值问题的子程序，对纳齐谢姆-斯威格特方法来说还需求解各边界函数对初值的一阶导数，其中也要用到求解常微分方程的子程序。这里选用Matlab函数ode45()来求解常微分方程。其调用格式为: [t,y]=ode45(@func,[tspan],y0,options,p1,p2…) 对本问题来说，用函数ode45()求解式(12-15)、式(12-16)时，先将该两式分解为5个一阶常微分方程，对布拉修斯方程其形式为: ,yy,,12,,yy,23,,,yyy,,,0.5 (12-33) ,313 ,,yy,45, ,,ypryy,,,-0.5515, 各边界函数对初值的一阶导数用向前差分表示。远端边界值及其一阶偏导数用矩阵形 9 式表示: ,,,fff()()(),,,，，xy,,,,,,,,fff()()(),,,xy,, (12-34) fn,,,,,,()()(),,,11,1,xy,, ,,,,,()()(),,,,,,11,1,xy，， 计算式(12-34)的函数为fn=fun_d_blasiusut(x,Pr,xspan,h)，其中x为假设的初值;Pr 为流体的普朗特数;xspan为无量纲离壁距离或相似变量取值区间，h为初值变化步长。对 m,0的楔形流动、自然对流等层流边界层流动与传热相似变换微分方程的求解只需重新 定义式(12-28)，并对远端边界条件做相应调整，也很容易用下面的程序求解。 用纳齐谢姆-斯威格特方法求解布拉修斯方程的完整Matlab程序如下: function blasius_N_S clc; clear all; close all; x0=[0.1,0.3]; etamax=8; pr=0.7;xspan=[0 etamax];h=0.01; eps=1e-9; n=1; xx1=x0+1; while(norm(xx1-x0)>eps)&(n<=100) f=fun_d_blasiusut(x0,pr,xspan,h); x1=f(1,2);x2=f(2,2); x3=f(3,2);x4=f(4,2); y1=f(1,3);y2=f(2,3); y3=f(3,3);y4=f(4,3); xx=x1*x1+x2*x2+x3*x3+x4*x4; yy=y1*y1+y2*y2+y3*y3+y4*y4; xy=x1*y1+x2*y2+x3*y3+x4*y4; cx=(1-f(1,1))*x1-f(2,1)*x2+(1-f(3,1))*x3-f(4,1)*x4; cy=(1-f(1,1))*y1-f(2,1)*y2+(1-f(3,1))*y3-f(4,1)*y4; qq=xx*yy-xy*xy; dx=(yy*cx-xy*cy)/qq;dy=(xx*cy-xy*cx)/qq; xx1=x0; x0(1)=x0(1)+dx;x0(2)=x0(2)+dy; n=n+1; 10 end y0=[0 0 x0(1) 0 x0(2)]; xspan=[0,etamax]; [x,y]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); figure(2) plot(x,y(:,2),'--r',x,y(:,4)) disp(double(y)); eer=abs(f(1,1)-1)+abs(f(2,1))+abs(f(3,1)-1)+abs(f(4,1)) err1=abs(f(1,1)-1) err3=abs(f(3,1)-1) err4=norm(xx1-x0) disp([n,xx1]); fprintf('%d %d %d %d %d',y(1,:)) function fn=fun_d_blasiusut(x,pr,xspan,h) x1=x; dx=h; dy=h; y0=[0;0;x1(1);0;x1(2)]; [x,y]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); a1=y(end,2); b1=y(end,3); c1=y(end,4); d1=y(end,5); x2=x1+h; y0=[0;0;x2(1);0;x1(2)]; [x,y]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); a2=y(end,2); b2=y(end,3); c2=y(end,4); d2=y(end,5); x2=x1+h; y0=[0;0;x1(1);0;x2(2)]; [x,y]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); a3=y(end,2); b3=y(end,3); c3=y(end,4); d3=y(end,5); ax=(a2-a1)/dx;ay=(a3-a1)/dy; 11 bx=(b2-b1)/dx;by=(b3-b1)/dy; cx=(c2-c1)/dx;cy=(c3-c1)/dy; dx=(d2-d1)/dx;dy=(d3-d1)/dy; fn=[a1 ax ay b1 bx by c1 cx cy d1 dx dy]; fn=fn+1e-7; 5. 应用Matlab软件的fsolve函数求解 fsolve函数先按方程求根计算，若不满足收敛条件给出函数极小值。 %blasius_ut_fsolve.m clc; clear all; close all; pr=[0.07 0.7 7]; etamax=[1.5 8 8]; x0=[0.3,0.3]; for k=1:3 slopew=fsolve('blasiust2',x0,[],pr(k),etamax(k)); y0=[0 0 slopew(1),0 slopew(2)]; xspan=[0,etamax(k)]; [eta,ff]=ode45(@(x,y0)blasiust(x,y0,pr(k)),xspan,y0); plot(eta,1-ff(:,4),'-') hold on end display([eta,1-ff(:,4)]) axis([0 5 0 1]) xlabel('\eta') ylabel('1-\theta') grid on ff 6. 应用Matlab软件的fminsearch函数求解 Matlab fminsearch函数，该函数优化算法采用单纯形法。 %blasisu_ut_fminsearch.m clc; clear; pr=0.7; etamax=20; 12 xm=5; x0=[1 0.7]; figure(1); slopew=fminsearch('blasiust2',x0,[],pr,etamax); slopew y0=[0 0 slopew(1) 0 slopew(2)]; xspan=[0 etamax]; [eta ff]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); subplot(2,1,1); plot(eta,ff(:,2),'-k',eta,ff(:,4),'--k'); subplot(2,1,2); plot(eta,ff(:,3),'--k',eta,ff(:,5),'-k'); ff 7. 应用Matlab软件的fminunc函数求解 Matlab fminunc函数，该函数优化算法采用拟Newton法。优化计算中Newton法的基 本格式为: ,1(1)()()()，kkkk ，，xxGx,,,f()，， ()k()k式中:为Hessian矩阵，为目标函数。微分方程求解采用四阶Runge-KuttaGf()x 法。 %blasius_ut_fmincun.m clc; clear; close all; pr=0.7; etamax=20; xm=5; x0=[1 0.7]; figure(1); opt.TolX=1e-8; slopew=fminunc('blasiust2',x0,[],pr,etamax); slopew y0=[0 0 slopew(1) 0 slopew(2)]; xspan=[0 etamax]; [eta ff]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); subplot(2,1,1); plot(eta,ff(:,2),'-k',eta,ff(:,4),'--k'); 13 subplot(2,1,2); plot(eta,ff(:,3),'--k',eta,ff(:,5),'-k'); ff 8. 应用Matlab软件的ga函数求解 function blasius_ga.m clc; clear all; x=[1 1]; numberOfVariables = 2; pr=0.7; eta_max=8; [x,fval] =ga(@(x)blasiust2(x,pr,eta_max),numberOfVariables) x,fval y0=[0,0,x(1),0,x(2)]; xspan=[0,eta_max]; [eta,ff]=ode45(@(x,y0)blasiust(x,y0,pr),xspan,y0); ff 14