Homework5.docx

1、Homework5Homework5Q: Solution:a. Let F(a,b,c)= .In order to find a,b,c adapt to minimize F.We can find the derivatives of F.=-Thus, we can get the aa= =2=0=2c=0=0Put a into this equation, we can get=0So we can get the nonlinear equation-We can also get the third equation like this- b. For this equat

2、ion and the two nonlinear equations, we can define six variables to replace these complex polynomials 1、t1=, t2=, t3=,t4=, t5=, t6=So these equations can be written asa=t1/t2;0=t1*t3-t2*t4;0=t1*t5-t2*t6;As for this problem, the Jacobi Matrix is 2X2. Its difficult to find the derivative of these form

3、ulas in MATLAB, so we can code them up directly. For example, for t1, we use these codes:function t1,dbt1,dct1=fun1(b,c,x,y,w) t1=0;dbt1=0;dct1=0; n=size(x,2);for i=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;endt1 is the equation,

4、dbt1 is, dct1 is. So its convenient to calculate for the Jacobi Matrix .The others can be written as:function t2,dbt2,dct2=fun2(b,c,x) t2=0;dbt2=0;dct2=0; n=size(x,2);for i=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);endfunction t3,dbt3,dct3=fun3(

5、b,c,x) t3=0;dbt3=0;dct3=0; n=size(x,2);for i=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);endfunction t4,dbt4,dct4=fun4(b,c,x,y,w) t4=0;dbt4=0;dct4=0; n=size(x,2);for i=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

6、+2); dct4=dct4-log(x(i)-b)*w(i)*y(i)/(x(i)-b)(c+1);endfunction t5, dbt5, dct5=fun5(b,c,x) t5=0;dbt5=0;dct5=0; n=size(x,2);for i=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);endfunction t6, dbt6, dct6=fun6(b,c,x, y, w) t6=0

7、;dbt6=0;dct6=0; n=size(x,2);for i=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;endSo the Jacobi Matrix can be written as: J(1,1)=dbt1*t3+t1*dbt3-dbt4*t2-t4*dbt2; J(1,2)=dct1*t3+t1*dct3-dct4*t2-t4*dct2; J(2

8、,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;So the Jacobi Matrix codes are: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);

9、 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;Then

10、we can start solving our problem1 Newtons methodfunction a,b,c,k,accuracy=zqhnew(x,y,b,c,N,TOL)w=log(y);k=1;while kN 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); if norm(b1-b c1-c,2) y=2.4 3.8 4.75 21.60; x=31.8 31.5 31.2 30.2; b=30; c=2; N=100; TOL=10-5; zqhne

11、w(x,y,b,c,N,TOL)The procedure was successful after 7 iterationsb = 26.7702c = 8.4518accuracy = 9.8830e-07ans = 2.2179e+062 Quasi-Newtons methodfunction 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);while kN b1=X1(1,1);c1=X1(2,1); b0=X0(1,1);c

12、0=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-A1F1; accuracy=norm(S,2); if norm(S,2) y=2.4 3.8 4.75 21.60; x=31.8 31.5 31.2 30.2; b=26; c=8; N=100; TOL=10-5; zqhquasi(x,y,b,c,N,TOL)The procedure was successful after 23

13、iterationsb = 25.9848c = 10.1123accuracy = 6.2590e-06ans = 1.3812e+083 Method of Steep Decentfunction a,b,c,k,g1,accuracy=zqhsteep(x,y,b,c,N,TOL)w=log(y);k=1;while k=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; if alpha3TOL/2 disp(NO likel

14、y improvement); 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-

15、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); if norm(g-g1,inf) y=2.4 3.8 4.75 21.60;x=31.8 31.5 31.2 30.2;b=30;c=2;N=100;TOL=10-5; zqhsteep(x,y,b,c,N,TOL)The procedure was successful after 5 iterationsb = 27.025

16、9c = 3.8817accuracy = 6.3440e-08ans = 5.0855e+034 Broyden Methodfunction 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);while kN Q=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*

17、V; b2=b1+s(1,1);c2=c1+s(2,1); accuracy=norm(c2-c1,inf); if norm(b2-b1 c2-c1,2) y=2.4 3.8 4.75 21.60;x=31.8 31.5 31.2 30.2;b=30;c=2;N=100;TOL=10-5;zqhbro(x,y,b,c,N,TOL)The procedure was successful after 11 iterationsb = 25.9848c = 10.1123accuracy = 6.2579e-06ans = 1.3812e+085 MATLAB Built-in Function

18、s(1) fminuncfunction zqhmm2ini_BC=30;2;options=optimset(Display,iter,Tolfun,10-7);opt_BC=fminunc(f,ini_BC,options) function fmin=f(BC)x=31.8 31.5 31.2 30.2;y=2.4 3.8 4.75 21.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;for i=1:4 fmin=

19、fmin+(w(i)*y(i)-a/(x(i)-BC(1,1)BC(2,1)2;end zqhmm2警告: Gradient must be provided for trust-region algorithm; using line-search algorithm instead. In fminunc at 382 In zqhmm2 at 4 First-order Iteration Func-count f(x) Step-size optimality 0 3 959.142 774 1 6 87.3578 0.0012913 290 2 9 18.1283 1 106 3 1

20、2 2.7427 1 35.9 4 15 1.44396 1 13.4 5 18 1.24402 1 0.879 6 21 1.24319 1 0.332 7 24 1.24292 1 0.398 8 27 1.24066 1 0.765 9 30 1.23626 1 1.65 10 33 1.22349 1 3.15 11 36 1.19264 1 4.93 12 39 1.11908 1 6.24 13 42 0.989803 1 5.52 14 45 0.877007 1 3.64 15 48 0.820879 1 2.25 16 51 0.796293 1 1.35 17 54 0.7

21、87799 1 0.819 18 57 0.78564 1 0.537 19 60 0.7851 1 0.403 First-order Iteration Func-count f(x) Step-size optimality 20 63 0.78467 1 0.3 21 66 0.783874 1 0.138 22 69 0.78291 1 0.0394 23 72 0.782308 1 0.111 24 75 0.782129 1 0.0421 25 78 0.782114 1 0.000993 26 81 0.782114 1 4.51e-06 Local minimum found

22、.Optimization completed because the size of the gradient is less thanthe selected value of the function tolerance.opt_BC = 26.76988.4526(2)fminsearchfunction zqhmmini_BC=30;2;options=optimset(Display,iter,Tolfun,10-7);opt_BC=fminsearch(f,ini_BC,options) function fmin=f(BC)x=31.8 31.5 31.2 30.2;y=2.4

23、 3.8 4.75 21.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;for i=1:4 fmin=fmin+(w(i)*y(i)-a/(x(i)-BC(1,1)BC(2,1)2;end zqhmm Iteration Func-count min f(x) Procedure 0 1 959.142 1 3 75.9237 initial simplex 2 5 52.7606 reflect 3 7 52.7606

24、 contract inside 4 8 52.7606 reflect 5 10 15.4821 contract inside 6 12 7.21226 contract inside 7 14 7.21226 contract outside 8 16 2.54277 contract inside 9 17 2.54277 reflect 10 19 2.54277 contract inside 11 21 2.3962 contract inside 12 23 2.21803 contract inside 13 24 2.21803 reflect 14 26 2.20942

25、contract inside 15 28 2.19528 contract inside 16 30 2.19326 reflect 17 32 2.17331 reflect 18 33 2.17331 reflect 19 35 2.15207 reflect 20 37 2.15207 contract inside 21 39 2.1336 expand 22 41 2.09239 expand 23 42 2.09239 reflect 24 44 2.00384 expand 25 45 2.00384 reflect 26 47 1.92683 expand 27 49 1.78149 e