Homework5.docx
《Homework5.docx》由会员分享，可在线阅读，更多相关《Homework5.docx（21页珍藏版）》请在冰豆网上搜索。
Homework5
Homework5
Q:
Solution:
a.LetF（a,b,c）=
.
Inordertofinda,b,cadapttominimizeF.
WecanfindthederivativesofF.
=
=-
Thus,wecangetthea
a=
=2
=0
=2c
=0
=0
Putaintothisequation,wecanget
=0
Sowecangetthenonlinearequation
-
Wecanalsogetthethirdequationlikethis
-
b.Forthisequationandthetwononlinearequations,wecandefinesixvariablestoreplacethesecomplexpolynomials
1、t1=
t2=
t3=
t4=
t5=
t6=
Sotheseequationscanbewrittenas
a=t1/t2;
0=t1*t3-t2*t4;
0=t1*t5-t2*t6;
Asforthisproblem,theJacobiMatrixis2X2.It’sdifficulttofindthederivativeoftheseformulasinMATLAB,sowecancodethemupdirectly.Forexample,fort1,weusethesecodes:
function[t1,dbt1,dct1]=fun1（b,c,x,y,w）
t1=0;dbt1=0;dct1=0;
n=size（x,2）;
fori=1:
n
t1=t1+w（i）*y（i）/（x（i）-b）^c;
dbt1=dbt1+c*（w（i）*y（i））/（x（i）-b）^（c+1）;
dct1=dct1-（w（i）*y（i）*log（x（i）-b））/（x（i）-b）^c;
end
t1istheequation,dbt1is
dct1is
.Soit’sconvenienttocalculatefortheJacobiMatrix.
Theotherscanbewrittenas:
function[t2,dbt2,dct2]=fun2（b,c,x）
t2=0;dbt2=0;dct2=0;
n=size（x,2）;
fori=1:
n
t2=t2+1/（x（i）-b）^（2*c）;
dbt2=dbt2+2*c/（x（i）-b）^（2*c+1）;
dct2=dct2-（2*log（x（i）-b））/（x（i）-b）^（2*c）;
end
function[t3,dbt3,dct3]=fun3（b,c,x）
t3=0;dbt3=0;dct3=0;
n=size（x,2）;
fori=1:
n
t3=t3+1/（x（i）-b）^（2*c+1）;
dbt3=dbt3+（2*c+1）/（x（i）-b）^（2*c+2）;
dct3=dct3-2*log（x（i）-b）/（x（i）-b）^（2*c+1）;
end
function[t4,dbt4,dct4]=fun4（b,c,x,y,w）
t4=0;dbt4=0;dct4=0;
n=size（x,2）;
fori=1:
n
t4=t4+（w（i）*y（i））/（x（i）-b）^（c+1）;
dbt4=dbt4+（c+1）*w（i）*y（i）/（x（i）-b）^（c+2）;
dct4=dct4-log（x（i）-b）*w（i）*y（i）/（x（i）-b）^（c+1）;
end
function[t5,dbt5,dct5]=fun5（b,c,x）
t5=0;dbt5=0;dct5=0;
n=size（x,2）;
fori=1:
n
t5=t5+log（x（i）-b）/（x（i）-b）^（2*c）;
dbt5=dbt5+（2*c*log（x（i）-b）-1）/（x（i）-b）^（2*c+1）;
dct5=dct5-2*log（x（i）-b）^2/（x（i）-b）^（2*c）;
end
function[t6,dbt6,dct6]=fun6（b,c,x,y,w）
t6=0;dbt6=0;dct6=0;
n=size（x,2）;
fori=1:
n
t6=t6+（log（x（i）-b）*w（i）*y（i））/（x（i）-b）^c;
dbt6=dbt6+w（i）*y（i）*（c*log（x（i）-b）-1）/（x（i）-b）^（c+1）;
dct6=dct6-w（i）*y（i）*log（x（i）-b）^2/（x（i）-b）^c;
end
SotheJacobiMatrixcanbewrittenas:
J（1,1）=dbt1*t3+t1*dbt3-dbt4*t2-t4*dbt2;
J（1,2）=dct1*t3+t1*dct3-dct4*t2-t4*dct2;
J（2,1）=dbt1*t5+t1*dbt5-dbt6*t2-t6*dbt2;
J（2,2）=dct1*t5+t1*dct5-dct6*t2-t6*dct2;
And
F（1,1）=t1*t3-t4*t2;
F（2,1）=t1*t5-t6*t2;
SotheJacobiMatrixcodesare:
function[F,J]=realfun（b,c,x,y,w）
F=zeros（2,1）;
J=zeros（2,2）;
[t1,dbt1,dct1]=fun1（b,c,x,y,w）;
[t2,dbt2,dct2]=fun2（b,c,x）;
[t3,dbt3,dct3]=fun3（b,c,x）;
[t4,dbt4,dct4]=fun4（b,c,x,y,w）;
[t5,dbt5,dct5]=fun5（b,c,x）;
[t6,dbt6,dct6]=fun6（b,c,x,y,w）;
F（1,1）=t1*t3-t4*t2;
F（2,1）=t1*t5-t6*t2;
J（1,1）=dbt1*t3+t1*dbt3-dbt4*t2-t4*dbt2;
J（1,2）=dct1*t3+t1*dct3-dct4*t2-t4*dct2;
J（2,1）=dbt1*t5+t1*dbt5-dbt6*t2-t6*dbt2;
J（2,2）=dct1*t5+t1*dct5-dct6*t2-t6*dct2;
Thenwecanstartsolvingourproblem
1Newton’smethod
function[a,b,c,k,accuracy]=zqhnew（x,y,b,c,N,TOL）
w=log（y）;
k=1;
whilek[F,J]=realfun（b,c,x,y,w）;
h=J\（-F）;
b1=b+h（1,1）;
c1=c+h（2,1）;
accuracy=norm（c1-c,inf）;
ifnorm（[b1-bc1-c],2）fprintf（'Theprocedurewassuccessfulafter%diterations',k）;
b
c
accuracy
[t1,dbt1,dct1]=fun1（b,c,x,y,w）;
[t2,dbt2,dct2]=fun2（b,c,x）;
a=t1/t2;
return
else
b=b1;c=c1;
k=k+1;
end
end
fprintf（'Methodfailedafter%diterations',N）;
end
>>y=[2.43.84.7521.60];
>>x=[31.831.531.230.2];
>>b=30;
>>c=2;
>>N=100;
>>TOL=10^-5;
>>zqhnew（x,y,b,c,N,TOL）
Theprocedurewassuccessfulafter7iterations
b=
26.7702
c=
8.4518
accuracy=
9.8830e-07
ans=
2.2179e+06
2Quasi-Newton’smethod
function[a,b,c,k,accuracy]=zqhquasi（x,y,b,c,N,TOL）
w=log（y）;
k=1;
[F0,J0]=realfun（b,c,x,y,w）;
A0=J0;
X0=[b;c];
X1=X0-J0\（F0）;
whilekb1=X1（1,1）;c1=X1（2,1）;
b0=X0（1,1）;c0=X0（2,1）;
[F1,J1]=realfun（b1,c1,x,y,w）;
[F0,J0]=realfun（b0,c0,x,y,w）;
S=X1-X0;Y=F1-F0;
A1=A0+（Y-A0*S）*S'/（（S（1,1））^2+（S（2,1））^2）;
X2=X1-A1\F1;
accuracy=norm（S,2）;
ifnorm（S,2）fprintf（'Theprocedurewassuccessfulafter%diterations',k）;
b=X1（1,1）
c=X1（2,1）
accuracy
[t1,dbt1,dct1]=fun1（b,c,x,y,w）;
[t2,dbt2,dct2]=fun2（b,c,x）;
a=t1/t2;
return
else
X0=X1;A0=A1;X1=X2;
k=k+1;
end
end
fprintf（'Methodfailedafter%diterations',N）;
end
>>y=[2.43.84.7521.60];
>>x=[31.831.531.230.2];
>>b=26;
>>c=8;
>>N=100;
>>TOL=10^-5;
>>zqhquasi（x,y,b,c,N,TOL）
Theprocedurewassuccessfulafter23iterations
b=
25.9848
c=
10.1123
accuracy=
6.2590e-06
ans=
1.3812e+08
3MethodofSteepDecent
function[a,b,c,k,g1,accuracy]=zqhsteep（x,y,b,c,N,TOL）
w=log（y）;
k=1;
whilek[F,J]=realfun（b,c,x,y,w）;
g1=F（1,1）^2+F（2,1）^2;
Z=2*J'*F;
z0=norm（Z,2）;
ifz0==0
disp（'Zerogradient'）;
k;b;c;g1;
[t1,dbt1,dct1]=fun1（b,c,x,y,w）;
[t2,dbt2,dct2]=fun2（b,c,x）;
a=t1/t2;
break;
end
Z=Z/z0;
alpha1=0;
alpha3=1;
b1=b-alpha3*Z（1,1）;
c1=c-alpha3*Z（2,1）;
[F,J]=realfun（b1,c1,x,y,w）;
g3=F（1,1）^2+F（2,1）^2;
whileg3>=g1
alpha3=alpha3/2;
b3=b-alpha3*Z（1,1）;
c3=c-alpha3*Z（2,1）;
[F,J]=realfun（b3,c3,x,y,w）;
g3=F（1,1）^2+F（2,1）^2;
ifalpha3disp（'NOlikelyimprovement'）;
k;b;c;g1;
[t1,dbt1,dct1]=fun1（b,c,x,y,w）;
[t2,dbt2,dct2]=fun2（b,c,x）;
a=t1/t2;
break;
end
end
alpha2=alpha3/2;
b2=b-alpha2*Z（1,1）;
c2=c-alpha2*Z（2,1）;
[F,J]=realfun（b2,c2,x,y,w）;
g2=F（1,1）^2+F（2,1）^2;
h1=（g2-g1）/alpha2;
h2=（g3-g2）/（alpha3-alpha2）;
h3=（h2-h1）/alpha3;
alpha0=0.5*（alpha2-h1/h3）;
b0=b-alpha0*Z（1,1）;
c0=c-alpha0*Z（2,1）;
[F,J]=realfun（b0,c0,x,y,w）;
g0=F（1,1）^2+F（2,1）^2;
g=g0;
bx=b0;cx=c0;
accuracy=norm（g-g1）;
ifnorm（g-g1,inf）fprintf（'Theprocedurewassuccessfulafter%diterations',k）;
b=bx
c=cx
accuracy
[t1,dbt1,dct1]=fun1（b,c,x,y,w）;
[t2,dbt2,dct2]=fun2（b,c,x）;
a=t1/t2;
break;
end
b=bx;
c=cx;
k=k+1;
end
>>y=[2.43.84.7521.60];
x=[31.831.531.230.2];
b=30;
c=2;
N=100;
TOL=10^-5;
>>zqhsteep（x,y,b,c,N,TOL）
Theprocedurewassuccessfulafter5iterations
b=
27.0259
c=
3.8817
accuracy=
6.3440e-08
ans=
5.0855e+03
4BroydenMethod
function[a,b,c,k,accuracy]=zqhbro（x,y,b,c,N,TOL）
w=log（y）;
k=2;
[F0,J0]=realfun（b,c,x,y,w）;
A0=J0;V=F0;
A=inv（A0）;
s=-A*V;
b1=b+s（1,1）;c1=c+s（2,1）;
whilekQ=V;
[F1,J1]=realfun（b1,c1,x,y,w）;
V=F1;P=V-Q;
z=-A*P;
p=-s'*z;
u=s'*A;
A=A+1/p*（s+z）*u;
s=-A*V;
b2=b1+s（1,1）;c2=c1+s（2,1）;
accuracy=norm（c2-c1,inf）;
ifnorm（[b2-b1c2-c1],2）fprintf（'Theprocedurewassuccessfulafter%diterations',k）;
b=b1
c=c1
accuracy
[t1,dbt1,dct1]=fun1（b,c,x,y,w）;
[t2,dbt2,dct2]=fun2（b,c,x）;
a=t1/t2;
return
else
b1=b2;c1=c2;
k=k+1;
end
end
fprintf（'Methodfailedafter%diterations',N）;
end
>>y=[2.43.84.7521.60];
x=[31.831.531.230.2];
b=30;
c=2;
N=100;
TOL=10^-5;
zqhbro（x,y,b,c,N,TOL）
Theprocedurewassuccessfulafter11iterations
b=
25.9848
c=
10.1123
accuracy=
6.2579e-06
ans=
1.3812e+08
5MATLABBuilt-inFunctions
（1）fminunc
functionzqhmm2
ini_BC=[30;2];
options=optimset（'Display','iter','Tolfun',10^-7）;
[opt_BC]=fminunc（@f,ini_BC,options）
function[fmin]=f（BC）
x=[31.831.531.230.2];
y=[2.43.84.7521.60];
w=log（y）;
[t1,t1db,t1dc]=fun1（BC（1,1）,BC（2,1）,x,y,w）;
[t2,t2db,t2dc]=fun2（BC（1,1）,BC（2,1）,x）;
a=t1/t2;
fmin=0;
fori=1:
4
fmin=fmin+（w（i）*y（i）-a/（x（i）-BC（1,1））^BC（2,1））^2;
end
>>zqhmm2
警告:
Gradientmustbeprovidedfortrust-regionalgorithm;
usingline-searchalgorithminstead.
>Infminuncat382
Inzqhmm2at4
First-order
IterationFunc-countf（x）Step-sizeoptimality
03959.142774
1687.35780.0012913290
2918.12831106
3122.7427135.9
4151.44396113.4
5181.2440210.879
6211.2431910.332
7241.2429210.398
8271.2406610.765
9301.2362611.65
10331.2234913.15
11361.1926414.93
12391.1190816.24
13420.98980315.52
14450.87700713.64
15480.82087912.25
16510.79629311.35
17540.78779910.819
18570.7856410.537
19600.785110.403
First-order
IterationFunc-countf（x）Step-sizeoptimality
20630.7846710.3
21660.78387410.138
22690.7829110.0394
23720.78230810.111
24750.78212910.0421
25780.78211410.000993
26810.78211414.51e-06
Localminimumfound.
Optimizationcompletedbecausethesizeofthegradientislessthan
theselectedvalueofthefunctiontolerance.
opt_BC=
26.7698
8.4526
（2）fminsearch
functionzqhmm
ini_BC=[30;2];
options=optimset（'Display','iter','Tolfun',10^-7）;
[opt_BC]=fminsearch（@f,ini_BC,options）
function[fmin]=f（BC）
x=[31.831.531.230.2];
y=[2.43.84.7521.60];
w=log（y）;
[t1,t1db,t1dc]=fun1（BC（1,1）,BC（2,1）,x,y,w）;
[t2,t2db,t2dc]=fun2（BC（1,1）,BC（2,1）,x）;
a=t1/t2;
fmin=0;
fori=1:
4
fmin=fmin+（w（i）*y（i）-a/（x（i）-BC（1,1））^BC（2,1））^2;
end
>>zqhmm
IterationFunc-countminf（x）Procedure
01959.142
1375.9237initialsimplex
2552.7606reflect
3752.7606contractinside
4852.7606reflect
51015.4821contractinside
6127.21226contractinside
7147.21226contractoutside
8162.54277contractinside
9172.54277reflect
10192.54277contractinside
11212.3962contractinside
12232.21803contractinside
13242.21803reflect
14262.20942contractinside
15282.19528contractinside
16302.19326reflect
17322.17331reflect
18332.17331reflect
19352.15207reflect
20372.15207contractinside
21392.1336expand
22412.09239expand
23422.09239reflect
24442.00384expand
25452.00384reflect
26471.92683expand
27491.78149e
升级会员 