最速下降法无约束最优化

1、MATLAB程序设计实践课程考核实践一、编程实现以下科学计算法，并举一例应用之。（参考书籍精通MATLAB科学计算，王正林等著，电子工业出版社，2009年）“最速下降法无约束最优化”最速下降法：解： 算法说明：最速下降法是一种沿着N维目标函数的负梯度方向搜索最小值的方法。原理：由高等数学知识知道任一点的负梯度方向是函数值在该点下降最快的方向，那么利用负梯度作为极值搜索方向，达到搜寻区间最速下降的目的。而极值点导数性质，知道该点的梯度=0，故而其终止条件也就是梯度逼近于0，也就是当搜寻区间非常逼近极值点时，即：当f(a）0推出f(a），f(a）即为所求。该方法是一种局部极值搜寻方法。函数的负梯度

2、表示如下：-g(x）=-f(x)=- 搜索步长可调整，通常记为k（第k次迭代中的步长）。该算法利用一维的线性搜索方法，如二次逼近法，沿着负梯度方向不断搜索函数的较小值，从而找到最优解。方法特点（1）初始值可任选，每次迭代计算量小，存储量少，程序简短。即使从一个不好的初始点出发，开始的几步迭代，目标函数值下降很快，然后慢慢逼近局部极小点。（2）任意相邻两点的搜索方向是正交的，它的迭代路径胃绕道逼近极小点。当迭代点接近极小点时，步长变得很小，越走越慢。（3）全局收敛，线性收敛，易产生扭摆现象而造成早停。算法步骤：最速下降法的基本求解流程如下：第一步迭代次数初始化为k=0，求出初始点的函数值f=f（

3、）。第二步迭代次数加1，即k=k+1，用一维线性搜索方法确定沿负梯度方向-的步长，其中=ArgMinaf（）。第三步沿着负梯度方向寻找下一个接近最小值的点，其中步长为，得到下一点的坐标为：。第四步如果，且f（）f（），那么就认为为所求的最小值点，并结束循环；否则，就跳到步骤二。流程图：-g(x）=-f(x)给定，ek=0开始。：minf（）=结束=ArgMinaf（）k=k+1是否题目：最速下降法求解无约束最优化问题实例。采用最速下降法求如下函数的最小值问题：f（x，y）=x(x-5-y)+y(y-4)即用最速下降法求解函数的最小值问题。解：需先求出该函数的梯度函数。可知其梯度函数为：g(x)

4、=(2x-5-y,-x+2y-4)。源程序代码如下：Opt_Steepest.m文件 %用最速下降法求最优化解；function xo,fo=Opt_Steepest(f,grad,x0,TolX,TolFun,dist0,MaxIter)%f:函数名；%grad:梯度函数；%x0：搜索初始值；%TolX:最优值点间的误差阈值；%TolFun:函数的误差阈值；%dist0:初始步长；%MaxIter:最大的迭代次数；%xo：最优化点值；%fo:函数在点xo处的函数值。%判断输入的变量数，设定一些变量为默认值if nargin<7 MaxIter=100; %最大的迭代次数默认为100en

5、dif nargin<6 dist0=10; %初始步长默认为10end if nargin<5 TolFun=1e-8; %函数值误差为1e-8endif nargin<4 TolX=1e-6; %自变量距离误差endx=x0;fx0=feval(f,x0);fx=fx0;dist=dist0;kmax1=25; %线性搜索法确定步长的最大搜索次数warning=0;%迭代计算求最优解for k=1:MaxIter g=feval(grad,x); g=g/norm(g); %求点x处的梯度 %线性搜索方法确定步长 dist=dist*2; fx1=feval(f,x-di

6、st*2*g); for k1=1:kmax1 fx2=fx1; fx1=feval(f,x-dist*g); if fx0>fx1+TolFun && fx1<fx2-TolFun %fx0>fx1<fx2,den=4*fx1-2*fx0-2*fx2;num=den-fx0+fx2; %二次逼近法 dist=dist*num/den; x=x-dist*g;fx=feval(f,x); %确定下一点 break; else dist=dist/2; end end if k1>=kmax1 warning=warning+1; %无法确定最优步长

7、 else warning=0; end if warning>=2|(norm(x-x0)<TolX&&abs(fx-fx0)<TolFun) break; end x0=x; fx0=fx;endxo=x;fo=fx;if k=MaxIter fprintf('Just best in %d iteration',MaxIter);end Q1.m文件f1004=inline('x(1)*(x(1)-5-x(2)+x(2)*(x(2)-4)','x'); %目标函数grad=inline('2*x(1

8、)-5-x(2),-x(1)+2*x(2)-4','x'); %目标函数的梯度函数x0=1 4;TolX=1e-4;TolFun=1e-9;MaxIter=100;dist0=1;xo,fo=Opt_Steepest(f1004,grad,x0,TolX,TolFun,dist0,MaxIter) 运行结果如下：由计算结果可知，当x=4.6667,y=4.3333时，函数f(x,y)=x(x-5-y)+y(y-4)取得最小值-20.3333。二编程解决以下科学计算和工程实际问题。简支梁受左半均匀分布载荷q及右边L/4处集中力偶M0作用（如下图1-1），求其弯矩、转角和挠

9、度。设L=2m，q=1000N/m，M0=900N*m，E=200*109N/m2,I=2*10-6m4.图1-1解题思路：首先对简支梁进行受力分析，受力分析图（如下图1-2）所示：图1-2从材料力学的知识可知道，由弯矩求转角要经过一次不定积分，而由转角求挠度又要经过一次不定积分，通常这是很麻烦而且容易出错的，而在MATLAB中，可用cumsum函数或cumtrapz函数作近似的不定积分，只要x取得足够密，其结果将相当准确，而程序非常简单。第一步：计算支反力设支座a和b处的支反力分别为Na和Nb，则据Ma=0，Fy=0得到平衡方程为：Nb=(q*L2/8+M0)/LNa=q*L/2-Nb第二步

10、：建立弯矩方程以截面c，d为分界面，将梁划分为ac，cd，db三段分别建立ac，cd，db三段对应的弯矩方程：M1=Na*x-q*x.2/2; 0xL/2 M2=Nb*(L-x)-M0; L/2x3L/4M3=Nb*(L-x); 3L/4xL第三步：建立挠曲轴近似微分方程并积分建立挠曲轴近似微分方程d2Y/dx2=M(x)/EI对M/EI积分，得转角A，再做一次积分，得挠度Y，每次积分都有一个待定的积分常数。A=(M)*dx/(E*I)+Ca= A0（x）+Ca， 此处设A0（x）= cumtrapz(M)*dx/(E*I) Y=(A)*dx+Cy=A0(x) *dx+Ca*x+Cy,此处设Y

11、0（x）=cumtrapz(A0)*dx第四步：确定相应的积分常数Ca，Cy由边界条件Y(0)=0，Y(L)=0确定Y(0)= Y0(0)+Cy=0Y(L)= Y0(L)+Ca*L+Cy=0即0 1；L 1 Ca；Cy=-Y0（0）；-Y0（L）Ca；Cy=0,1;L,1-Y0(0);-Y0（L);第五步：根据计算结果绘制弯矩、转角以及挠度图形源程序：L=2;q=1000;M0=900;E=200e9;I=2e-6; %输入已知参数Nb=(q*L2/8+M0)/L;Na=q*L/2-Nb; %求支反力x=linspace(0,L,101);dx=L/100; %linspace是线性空间向量M

12、1=Na*x(1:51)-q*x(1:51).2/2; %分三段用数组列出M数组M2=Nb*(L-x(52:76)-M0;M3=Nb*(L-x(77:101);M=M1,M2,M3; %写出完整的M数组A0=cumtrapz(M)*dx/(E*I); %由M积分求转角AY0=cumtrapz(A0)*dx; %由转角积分求挠度YC=0,1;L,1-Y0(1);-Y0(101); %由边界条件求积分常数Ca=C(1),Cy=C(2),A=A0+Ca;Y=Y0+Ca*x+Cy; %A为转角，Y为挠度，求出转角和挠度的完整数值subplot(3,1,1),plot(x,M),grid %绘制弯矩图形

13、subplot(3,1,2),plot(x,A),grid %绘制转角图形subplot(3,1,3),plot(x,Y),grid %绘制挠度图形 运行结果：开始输入已知参数L，q，M0，E，I求支反力Nb，Na使x=linspace(0,L,101)，dx=L/100，划分x为空间线性向量分三段用数组列出M数组，写出完整的M数组M=M1,M2,M3由M积分求转角A，由转角积分求挠度Y （用cumtrapz积分），由边界条件求积分常数求出转角和挠度的完整数值，A=A0+Ca;Y=Y0+Ca*x+Cy;分别绘制弯矩，转角，挠度的图形用subplot分别绘制弯矩，转角，挠度的图形并输出输出积分常

14、数Ca,Cy结束流程图：第三题流程图：开始输入已知的数据表作为样本；设置插值节点针对不同的方法选用相应的函数及格式将已知数据和插值节点代入求得插值节点处的函数值（1）A正弦值算法：x=0:pi/12:pi/2;y=0 0.2588 0.5000 0.7071 0.8660 0.9659 1.0000;xi=0:pi/180:pi/2;%三次样条差值yi=interp1(x,y,xi,'spline')%五次多项式拟合A=polyfit(x,y,5);yj=polyval(A,xi)运行结果：yi = Columns 1 through 11 0 0.0175 0.0349 0.

15、0524 0.0698 0.0872 0.1045 0.1219 0.1392 0.1564 0.1737 Columns 12 through 22 0.1908 0.2079 0.2249 0.2419 0.2588 0.2756 0.2923 0.3090 0.3255 0.3420 0.3583 Columns 23 through 33 0.3746 0.3907 0.4067 0.4226 0.4384 0.4540 0.4695 0.4848 0.5000 0.5150 0.5299 Columns 34 through 44 0.5446 0.5592 0.5736 0.587

16、8 0.6018 0.6157 0.6293 0.6428 0.6561 0.6691 0.6820 Columns 45 through 55 0.6947 0.7071 0.7193 0.7313 0.7431 0.7547 0.7660 0.7771 0.7880 0.7986 0.8090 Columns 56 through 66 0.8191 0.8290 0.8387 0.8480 0.8571 0.8660 0.8746 0.8829 0.8910 0.8987 0.9062 Columns 67 through 77 0.9135 0.9204 0.9271 0.9335 0

17、.9396 0.9454 0.9510 0.9563 0.9612 0.9659 0.9703 Columns 78 through 88 0.9744 0.9782 0.9817 0.9849 0.9878 0.9904 0.9927 0.9946 0.9963 0.9977 0.9987 Columns 89 through 91 0.9995 0.9999 1.0000yj = Columns 1 through 11 0.0000 0.0174 0.0349 0.0523 0.0697 0.0871 0.1045 0.1218 0.1391 0.1564 0.1736 Columns

18、12 through 22 0.1908 0.2079 0.2249 0.2419 0.2588 0.2756 0.2924 0.3090 0.3256 0.3420 0.3584 Columns 23 through 33 0.3746 0.3907 0.4067 0.4226 0.4384 0.4540 0.4695 0.4848 0.5000 0.5150 0.5299 Columns 34 through 44 0.5446 0.5592 0.5736 0.5878 0.6018 0.6157 0.6293 0.6428 0.6561 0.6691 0.6820 Columns 45

19、through 55 0.6946 0.7071 0.7193 0.7313 0.7431 0.7547 0.7660 0.7771 0.7880 0.7986 0.8090 Columns 56 through 66 0.8191 0.8290 0.8386 0.8480 0.8571 0.8660 0.8746 0.8829 0.8910 0.8988 0.9063 Columns 67 through 77 0.9135 0.9205 0.9272 0.9336 0.9397 0.9455 0.9510 0.9563 0.9612 0.9659 0.9703 Columns 78 thr

20、ough 88 0.9743 0.9781 0.9816 0.9848 0.9877 0.9902 0.9925 0.9945 0.9962 0.9975 0.9986 Columns 89 through 910.9994 0.9998 1.0000通过比较，两种方法得到的结果近似相等。B正切值算法：x=0:pi/12:5*pi/12;y=0 0.2679 0.5774 1.0000 1.7320 3.7320;xi=0:pi/180:5*pi/12;%三次样条差值yi=interp1(x,y,xi,'spline')%五次多项式拟合A=polyfit(x,y,5);yj=p

21、olyval(A,xi)运行结果：yi = Columns 1 through 11 0 0.0184 0.0365 0.0545 0.0724 0.0902 0.1079 0.1255 0.1431 0.1607 0.1784 Columns 12 through 22 0.1961 0.2138 0.2317 0.2497 0.2679 0.2863 0.3048 0.3236 0.3427 0.3620 0.3817 Columns 23 through 33 0.4017 0.4221 0.4429 0.4641 0.4858 0.5079 0.5305 0.5537 0.5774 0

22、.6017 0.6266 Columns 34 through 44 0.6520 0.6780 0.7046 0.7317 0.7593 0.7876 0.8163 0.8456 0.8754 0.9058 0.9367 Columns 45 through 55 0.9681 1.0000 1.0325 1.0658 1.1003 1.1364 1.1743 1.2145 1.2572 1.3028 1.3516 Columns 56 through 66 1.4041 1.4604 1.5211 1.5863 1.6565 1.7320 1.8131 1.9002 1.9936 2.09

23、37 2.2008 Columns 67 through 76 2.3152 2.4374 2.5675 2.7060 2.8532 3.0095 3.1752 3.3506 3.5361 3.7320yj = Columns 1 through 11 -0.0000 0.0235 0.0454 0.0658 0.0850 0.1032 0.1206 0.1375 0.1540 0.1701 0.1862 Columns 12 through 22 0.2022 0.2183 0.2345 0.2511 0.2679 0.2851 0.3028 0.3208 0.3394 0.3585 0.3

24、781 Columns 23 through 33 0.3982 0.4188 0.4400 0.4616 0.4838 0.5065 0.5297 0.5533 0.5774 0.6020 0.6270 Columns 34 through 44 0.6524 0.6783 0.7047 0.7315 0.7588 0.7867 0.8150 0.8440 0.8736 0.9039 0.9351 Columns 45 through 55 0.9670 1.0000 1.0341 1.0693 1.1060 1.1442 1.1841 1.2259 1.2699 1.3162 1.3652

25、 Columns 56 through 66 1.4171 1.4723 1.5310 1.5935 1.6604 1.7320 1.8087 1.8910 1.9793 2.0742 2.1762 Columns 67 through 762.2860 2.4040 2.5310 2.6677 2.8147 2.9727 3.1427 3.3253 3.5214 3.7320通过比较知，角度较小时五次多项式算得的值较大，角度增大则两种方法得到的结果近似相等。(2) %三次多项式插值x=1 4 9 16 25 36 49 64 81 100;y=1 2 3 4 5 6 7 8 9 10;xi=1:100;f=interp1(x,y,xi,'cubic')运行结果：f = Columns 1 through 11 1.0000 1.3729 1.7125 2.0000 2.2405 2.4551