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

《MATLAB程序设计实践》课程考核

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

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

算法步骤：最速下降法的基本求解流程如下：
第一步
迭代次数初始化为k=0，求出初始点的函数值f=f（）。
第二步
迭代次数加1，即k=k+1，用一维线性搜索方法确定沿负梯度方向-的步长，其中=ArgMinaf（）。
第三步
沿着负梯度方向寻找下一个接近最小值的点，其中步长为，得到下一点的坐标为：。
第四步
如果≈，且f（）≈f（），那么就认为为所求的最小值点，并结束循环；否则，就跳到步骤二。

流程图：
-g(x）=-▽f(x)
给定，e
k=0
开始
。
：minf（）

=
结束
=ArgMinaf（）
k=k+1
是
否

题目：
最速下降法求解无约束最优化问题实例。采用最速下降法求如下函数的最小值问题：
f（x，y）=x(x-5-y)+y(y-4)
即用最速下降法求解函数的最小值问题。
解：需先求出该函数的梯度函数。可知其梯度函数为：g(x)=(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; %最大的迭代次数默认为100
end
if nargin<6
dist0=10; %初始步长默认为10
end
if nargin<5
TolFun=1e-8; %函数值误差为1e-8
end
if nargin<4
TolX=1e-6; %自变量距离误差
end
x=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-dist*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; %无法确定最优步长
else
warning=0;
end
if warning>=2||(norm(x-x0)<TolX&&abs(fx-fx0)<TolFun)
break;
end
x0=x;
fx0=fx;
end
xo=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)-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），求其弯矩、转角与挠度。设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*L^2/8+M0)/L
Na=q*L/2-Nb
第二步：建立弯矩方程
以截面c，d为分界面，将梁划分为ac，cd，db三段
分别建立ac，cd，db三段对应的弯矩方程：
M1=Na*x-q*x.^2/2; 0≦x≦L/2
M2=Nb*(L-x)-M0; L/2≦x≦3L/4
M3=Nb*(L-x); 3L/4≦x≦L
第三步：建立挠曲轴近似微分方程并积分
建立挠曲轴近似微分方程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,此处设Y0（x）=cumtrapz(A0)*dx
第四步：确定相应的积分常数
Ca，Cy由边界条件Y(0)=0，Y(L)=0确定
Y(0)= Y0(0)+Cy=0
Y(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*L^2/8+M0)/L;Na=q*L/2-Nb; %求支反力
x=linspace(0,L,101);dx=L/100; %linspace是线性空间向量
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积分求转角A
Y0=cumtrapz(A0)*dx; %由转角积分求挠度Y
C=[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 %绘制弯矩图形
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分别绘制弯矩，转角，挠度的图形并输出
输出积分常数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.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.5878 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.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.0000

yj =

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 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 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 through 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 91

0.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=polyval(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.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.0937 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.7320

yj =

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.3781

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

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 76

2.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 2.6494 2.8292 3.0000 3.1636 3.3186

Columns 12 through 22

3.4661 3.6069 3.7422 3.8729 4.0000 4.1237 4.2435 4.3599 4.4730 4.5832 4.6907

Columns 23 through 33

4.7958 4.8988 5.0000 5.0993 5.1966 5.2921 5.3857 5.4777 5.5681 5.6570 5.7446

Columns 34 through 44

5.8309 5.9160 6.0000 6.0829 6.1647 6.2454 6.3249 6.4035 6.4810 6.5577 6.6334

Columns 45 through 55

6.7082 6.7823 6.8556 6.9281 7.0000 7.0712 7.1416 7.2113 7.2804 7.3487 7.4164

Columns 56 through 66

7.4835 7.5500 7.6159 7.6812 7.7459 7.8102 7.8739 7.9372 8.0000 8.0623 8.1242

Columns 67 through 77

8.1855 8.2464 8.3068 8.3668 8.4263 8.4854 8.5441 8.6024 8.6603 8.7178 8.7749

Columns 78 through 88

8.8317 8.8881 8.9442 9.0000 9.0555 9.1107 9.1655 9.2201 9.2744 9.3284 9.3821

Columns 89 through 99

9.4354 9.4884 9.5412 9.5935 9.6456 9.6973 9.7486 9.7996 9.8502 9.9005 9.9505

Column 100

10.0000

27 / 27