MATLAB程序设计实践课程考核.doc

《MATLAB程序设计实践课程考核.doc》由会员分享，可在线阅读，更多相关《MATLAB程序设计实践课程考核.doc（15页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

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

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

3、，用一维线性搜索方法确定沿负梯度方向-步长，其中=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)=(2x-5-y,-x+2y-4)。源程序代码如下

4、：文件 %用最速下降法求最优化解；function xo,fo=Opt_Steepest(f,grad,x0,TolX,TolFun,dist0,MaxIter)%f:函数名；%grad:梯度函数；%x0：搜索初始值；%TolX:最优值点间误差阈值；%TolFun:函数误差阈值；%dist0:初始步长；%MaxIter:最大迭代次数；%xo：最优化点值；%fo:函数在点xo处函数值。%判断输入变量数，设定一些变量为默认值if nargin7 MaxIter=100; %最大迭代次数默认为100endif nargin6 dist0=10; %初始步长默认为10end if nargin5 To

5、lFun=1e-8; %函数值误差为1e-8endif narginfx1+TolFun & fx1fx1=kmax1 warning=warning+1; %无法确定最优步长 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);endQ1.m文件f1004=inline(x(1)*(x(1)-5-x(2)+x(2)*(x(2

6、)-4),x); %目标函数grad=inline(2*x(1)-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作用（如下

7、图1-1），求其弯矩、转角与挠度。设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=

8、q*L/2-Nb第二步：建立弯矩方程以截面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*

9、x+Cy,此处设Y0（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

10、是线性空间向量M1=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 %

11、绘制弯矩图形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分别绘制弯矩，转角，挠度的图形并输

12、出输出积分常数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 Columns 12 through 22

13、Columns 23 through 33 Columns 34 through 44 Columns 45 through 55 Columns 56 through 66 Columns 67 through 77 Columns 78 through 88 Columns 89 through 91yj = Columns 1 through 11 Columns 12 through 22 Columns 23 through 33 Columns 34 through 44 Columns 45 through 55 Columns 56 through 66 Columns 67

14、through 77 Columns 78 through 88 Columns 89 through 91通过比较，两种方法得到结果近似相等。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=polyval(A,xi)运行结果：yi = Columns 1 through 11 Columns 12 through 22 Columns 23 through

15、 33 Columns 34 through 44 Columns 45 through 55 Columns 56 through 66 Columns 67 through 76yj = Columns 1 through 11 Columns 12 through 22 Columns 23 through 33 Columns 34 through 44 Columns 45 through 55 Columns 56 through 66 Columns 67 through 76通过比较知，角度较小时五次多项式算得值较大，角度增大则两种方法得到结果近似相等。(2) %三次多项式插值

16、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 Columns 12 through 22 Columns 23 through 33 Columns 34 through 44 Columns 45 through 55 Columns 56 through 66 Columns 67 through 77 Columns 78 through 88 Columns 89 through 99 Column 100第 15 页