对于本次编程编程作业，小组运用matlab 和c++两种程序对平板边界层问题和绕过楔形体边界层流动问题进行分析研究。

以下是运用matlab 解决问题的过程。

一、 平板边界层问题该问题可以归结为在已知边界层条件下解一个高阶微分方程,即解0''5.0'''=+ff f 。

Matlab 提供了解微分方程的方法，运用换元法将高阶微分方程降阶，然后运用“ode45”函数进行求解。

函数其难点在于如何将边界条件中1',→∞→f η运用好，由四阶龙格-库塔方法知其核心是换元试算匹配，故在运用函数时通过二分法实现1',→∞→f η是可行的。

程序如下：第一问m 函数function dy = rigid(t,y) dy = zeros(3,1); dy(1) = y(2); dy(2) = y(3);dy(3) = -0.5*y(1)*y(3);%第一问main 程序[T,Y] = ode45('rigid',[0 5],[0 0 0]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') %二分法试算f ’’的初始值以满足f ’趋向无穷时的边界条件，图像上可以清晰看出f ’无穷时的结果 >> [T,Y] = ode45('rigid',[0 5],[0 0 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> [T,Y] = ode45('rigid',[0 5],[0 0 0.5]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> [T,Y] = ode45('rigid',[0 5],[0 0 0.25]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')>> [T,Y] = ode45('rigid',[0 5],[0 0 0.375]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> grid on>> [T,Y] = ode45('rigid',[0 5],[0 0 0.3125]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> grid on>> [T,Y] = ode45('rigid',[0 5],[0 0 0.34375]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on>> [T,Y] = ode45('rigid',[0 5],[0 0 0.328125]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on>> [T,Y] = ode45('rigid',[0 10],[0 0 0.328125]);%当f ’’为0.328125时，逼近结果已经很好，在0到5的变化范围内已经非常接近精确解 plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on>> [T,Y] = ode45('rigid',[0 5],[0 0 0.335975]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on选取f’’=0.335975时的数据展示，T代表η的变化，Y的第一二三列代表ff,f,'''为形象展示所得结果，如图-1所示：二、楔形绕流边界层问题解决绕过楔形的边界层流动问题与平板问题的不同之处在于微分方程变为0α，其中α可取1，并不失其普遍性，由于β+ffβ+ff'''2=)'1(''-小于-0.199时会发生分离，该微分方程不再适用，故β取值大于-0.199。

其解法与平板绕流类似，在平板绕流基础上增加β的变化即可。

程序如下：第二问m函数function dz = rigid1(t,z)dz = zeros(3,1);t=-0.199; %其中t表示β的变化，可由-0.199至无穷大dz(1) = z(2);dz(2) = z(3);dz(3) = -z(1)*z(3)-t*(1-z(2)^2);第二问main函数[T,Y] = ode45('rigid1',[0 5],[0 0 0]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') %二分法试算f’’的初始值以满足f’趋向无穷时的边界条件，图像上可以清晰看出f’无穷时的结果>> [T,Y] = ode45('rigid1',[0 5],[0 0 1]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')>> [T,Y] = ode45('rigid1',[0 5],[0 0 0.5]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')>> [T,Y] = ode45('rigid1',[0 5],[0 0 0.25]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')>> [T,Y] = ode45('rigid1',[0 5],[0 0 0.15]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')>> grid on>> [T,Y] = ode45('rigid1',[0 5],[0 0 0.05]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')>> grid on>> [T,Y] = ode45('rigid1',[0 5],[0 0 0.025]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')grid on%当f’’为0.025时，逼近结果已经很好，在0到5的变化范围内已经非常接近精确解plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')grid on>> [T,Y] = ode45('rigid1',[0 5],[0 0 0]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')grid on选取β=-0.199,f’’=0.025时的数据进行展示，程序代码中T代表η的变化，y(i)的代表''f的值。

第二个问题的数据如下表：f,f,'当β变化时，只需要更改主程序中t的值即可，再运用main程序解微分方程即可得到相应的解。

三、优化楔形绕流问题的算法初值的选取对于减少计算量有很大帮助，首先考虑运用拟合得到经验性公式，尝试如下：在β变化时按此式（其中t代表β）f''=-0.2877t^2+1.1892t+0.3915变化f’’相应的初值利于快速得到结果。

该式子是由二阶拟合多组数据得到的结果。

具体程序为x = [-0.1988 -0.19 -0.18 -0.16 -0.14 -0.10 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.80 1.00 1.20 1.60 2.00];y = [0 0.0860 0.1285 0.1905 0.2395 0.3191 0.4696 0.5870 0.6869 0.77680.8542 0.9277 0.9960 1.120 1.233 1.336 1.521 1.687];A = polyfit(x,y,2)z = polyval(A,x);plot(x,y,'k+',x,z,'r')%f''=0.2133*t^3-0.8160*t^2+1.4229*t+0.4189 三阶公式%f''=0.2826*t^5-1.4303*t^4+2.6049*t^3-2.1971*t^2+1.5037*t+0.4723 五阶公式图像如图-3因此，在β发生变化时只需要在主程序中相应的变化t值即可。

例如取β=1.2时，由f''=-0.2877t^2+1.1892t+0.3915式子可得初值为1.4030,更改m函数中t值为1.2，再调用main函数，并赋予f’’新的初值，如下：[T,Y] = ode45('rigid1',[0 5],[0 0 1.3199]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')grid on应用过程中发现f’并不能较快的趋近于1，即使运用四次拟合公式（f''=0.2826*t^5-1.4303*t^4+2.6049*t^3-2.1971*t^2+1.5037*t+0.4723）也不能较快的得到收敛的结果，其主要原因是精度达不到要求和拟合本身的误差过大。