最速下降法无约束最优化Word下载.doc
《最速下降法无约束最优化Word下载.doc》由会员分享,可在线阅读,更多相关《最速下降法无约束最优化Word下载.doc(17页珍藏版)》请在冰点文库上搜索。
否则,就跳到步骤二。
流程图:
-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处的函数值。
%%%%%%判断输入的变量数,设定一些变量为默认值
ifnargin<
7
MaxIter=100;
%最大的迭代次数默认为100
end
6
dist0=10;
%初始步长默认为10
end
5
TolFun=1e-8;
%函数值误差为1e-8
4
TolX=1e-6;
%自变量距离误差
x=x0;
fx0=feval(f,x0);
fx=fx0;
dist=dist0;
kmax1=25;
%线性搜索法确定步长的最大搜索次数
warning=0;
%%%%%迭代计算求最优解
fork=1:
MaxIter
g=feval(grad,x);
g=g/norm(g);
%求点x处的梯度
%%线性搜索方法确定步长
dist=dist*2;
fx1=feval(f,x-dist*2*g);
fork1=1:
kmax1
fx2=fx1;
fx1=feval(f,x-dist*g);
iffx0>
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
ifk1>
=kmax1
warning=warning+1;
%无法确定最优步长
else
warning=0;
ifwarning>
=2||(norm(x-x0)<
TolX&
abs(fx-fx0)<
TolFun)
break;
x0=x;
fx0=fx;
xo=x;
fo=fx;
ifk==MaxIter
fprintf('
Justbestin%diteration'
MaxIter);
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]'
%目标函数的梯度函数
x0=[14];
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
即[01;
L1][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;
分别绘制弯矩,转角,挠度的图形
用subplot分别绘制弯矩,转角,挠度的图形并输出
输出积分常数Ca,Cy
④流程图:
第三题
输入已知的数据表作为样本;
设置插值节点
针对不同的方法选用相应的函数及格式
将已知数据和插值节点代入
求得插值节点处的函数值
(1)
A.正弦值算法:
x=0:
pi/12:
pi/2;
y=[00.25880.50000.70710.86600.96591.0000];
xi=0:
pi/180:
%三次样条差值
yi=interp1(x,y,xi,'
spline'
)
%五次多项式拟合
A=polyfit(x,y,5);
yj=polyval(A,xi)
运行结果:
yi=
Columns1through11
00.01750.03490.05240.06980.08720.10450.12190.13920.15640.1737
Columns12through22
0.19080.20790.22490.24190.25880.27560.29230.30900.32550.34200.3583
Columns23through33
0.37460.39070.40670.42260.43840.45400.46950.48480.50000.51500.5299
Columns34through44
0.54460.55920.57360.58780.60180.61570.62930.64280.65610.66910.6820
Columns45through55
0.69470.70710.71930.73130.74310.75470.76600.77710.78800.79860.8090
Columns56through66
0.81910.82900.83870.84800.85710.86600.87460.88290.89100.89870.9062
Columns67through77
0.91350.92040.92710.93350.93960.94540.95100.95630.96120.96590.9703
Columns78through88
0.97440.97820.98170.98490.98780.99040.99270.99460.99630.99770.9987
Columns89through91
0.99950.99991.0000
yj=
0.00000.01740.03490.05230.06970.08710.10450.12180.13910.15640.1736
0.19080.20790.22490.24190.25880.27560.29240.30900.32560.34200.3584
0.69460.70710.71930.73130.74310.75470.76600.77710.78800.79860.8090
0.81910.82900.83860.84800.85710.86600.87460.88290.89100.89880.9063
0.91350.92050.92720.93360.93970.94550.95100.95630.96120.96590.9703
0.97430.97810.98160.98480.98770.99020.99250.99450.99620.99750.9986
0.99940.99981.0000
通过比较,两种方法得到的结果近似相等。
B.正切值算法:
5*pi/12;
y=[00.26790.57741.00001.73203.7320];
00.01840.03650.05450.07240.09020.10790.12550.14310.16070.1784
0.19610.21380.23170.24970.26790.28630.30480.32360.34270.36200.3817
0.40170.42210.44290.46410.48580.50790.53050.55370.57740.60170.6266
0.65200.67800.70460.73170.75930.78760.81630.84560.87540.90580.9367
0.96811.00001.03251.06581.10031.13641.17431.21451.25721.30281.3516
1.40411.46041.52111.58631.65651.73201.81311.90021.99362.09372.2008
Columns67through76
2.31522.43742.56752.70602.85323.00953.17523.35063.53613.7320
-0.00000.02350.04540.06580.08500.10320.12060.13750.15400.17010.1862
0.20220.21830.23450.25110.26790.28510.30280.32080.33940.35850.3781
0.39820.41880.44000.46160.48380.50650.52970.55330.57740.60200.6270
0.65240.67830.70470.73150.75880.78670.81500.84400.87360.90390.9351
0.96701.00001.03411.06931.10601.14421.18411.22591.26991.31621.3652
1.41711.47231.53101.59351.66041.73201.80871.89101.97932.07422.1762
2.28602.40402.53102.66772.81472.97273.14273.32533.52143.7320
通过比较知,角度较小时五次多项式算得的值较大,角度增大则两种方法得到的结果近似相等。
(2)
%三次多项式插值
x=[149162536496481100];
y=[12345678910];
xi=1:
100;
f=interp1(x,y,xi,'
cubic'
f=
1.00001.37291.71252.00002.24052.45512