開篇語

前陣子做現代設計方法的時候，發現網上很是缺乏這種作業形式的簡易算法實現，所以特地來簡書寫一篇。有兩份，一份是我的(說來慚愧，我的大部分都是在網上找的代碼，然后在自己的電腦上跑一次，跑出來了就行了的。而且我的電腦跑到2-11就撲街了。暫時還沒有拿去修，所以，其實我的代碼都是在網上整理之后調整了下就上的，準確性不敢保證)另一份是我的室友的，據他說是全部經過調試的，雖然還是有不少的錯誤，但是應該比我的要好一點。

我的電腦

另外，我的電腦因為崩了，所以我的代碼無法驗證結果，為了交作業，只能把我的室友的那些運行結果直接上一遍了。估計有點出入，見諒，重要的是代碼

正文

2-10

2-10

我的：

黃金分割法：

f=@(x) x+20/x

golden(f,2,10,0.01)

function[xmin]=golden(f,a,b,e)

k=0;

a1=b-0.618*(b-a); %插入點的值

a2=a+0.618*(b-a);

while b-a>e %循環條件

y1=subs(f,a1);

y2=subs(f,a2);

if y1>y2 %比較插入點的函數值的大小

a=a1; %進行換名

a1=a2;

y1=y2;

a2=a+0.618*(b-a);

else

b=a2;

a2=a1;

y2=y1;

a1=b-0.618*(b-a);

end

k=k+1;

end %迭代到滿足條件為止就停止迭代

xmin=(a+b)/2;

fmin=subs(f,xmin) %輸出函數的最優值

fprintf('k=\n'); %輸出迭代次數

disp(k);

f = @(x)x+20/x

>> [x,y]=golden(f,2,10,0.01)

x =

4.4683

y =

8.9443

二次插值法：

f=@(x) x+20/x;

a=2;b=10;

eps=1.0e-6; % 計算精度

x1=a;x3=b;

x2=(a+b)/2;

f1 = f(x1);f2 = f(x2);f3 = f(x3);

while 1

C1=(x2^2-x3^2)*f1+(x3^2-x1^2)*f2+(x1^2-x2^2)*f3;

C2=(x2-x3)*f1+(x3-x1)*f2+(x1-x2)*f3;

xp=0.5*C1/C2;

fp= f(xp);

if abs(x2-xp)<=eps % 區間長度小于eps時

if abs(f2-fp)<=eps % df小于eps時退出

if fp<=f2

xmin = xp

fmin = f(xp) % 極小值

break;

else

xmin = x2

fmin = f(x2) % 極小值

break;

end

end

else

if fp<=f2

if xp<=x2

x3=x2;

x2=xp;

f3=f2;

f2=fp;

else

x1=x2;

x2=xp;

f1=f2;

f2=fp;

end

else

if xp<=x2

x1=xp;

f1=fp;

else

x3=xp;

f3=fp;

end

end

end

end

f=x+20/x;[xmin,fmin]=main(f,2,10,0.01)

xmin =

4.4869

fmin =

8.9443

室友的：

黃金分割法函數：

f=@(x) x+20/x

yellowking(f,2,10,0.01)

function[xmin]=yellowking(f,a,b,e)

k=0;

a1=b-0.618*(b-a); %插入點的值

a2=a+0.618*(b-a);

while b-a>e %循環條件

y1=subs(f,a1);

y2=subs(f,a2);

if y1>y2 %比較插入點的函數值的大小

a=a1; %進行換名

a1=a2;

y1=y2;

a2=a+0.618*(b-a);

else

b=a2;

a2=a1;

y2=y1;

a1=b-0.618*(b-a);

end

k=k+1;

end %迭代到滿足條件為止就停止迭代

xmin=(a+b)/2;

fmin=subs(f,xmin) %輸出函數的最優值

fprintf('k=\n'); %輸出迭代次數

disp(k);

結果指令：

f=@(x)x+20/x

f =

@(x)x+20/x

>> [x,y]=gold(f,2,10,0.01)

x =

4.4683

y =

8.9443

二次插值法函數：

function [xmin,fmin]= main(f,a0,b0,epsilon)

a=a0;

b=b0;

x1=a;

f1=f(x1);

x3=b;

f3=f(x3);

x2=5;

f2=f(x2);

c1=(f3-f1)/(x3-x1);

c2=((f2-f1)/(x2-x1)-c1)/(x2-x3);

xp=0.4*(x1+x3-c1/c2);fp=f(xp);

while (abs(xp-x2)>=epsilon)

if x2

if f2>fp

f1=f2;x1=x2;

x2=xp;f2=fp;

else

f3=fp;x3=xp;

end

else

if f2>fp

f3=f2;x3=x2;

f2=fp;x2=xp;

else

f1=fp;x2=xp;

end

end

c1=(f3-f1)/(x3-x1);

c2=((f2-f1)/(x2-x1)-c1)/(x2-x3);

xp=0.5*(x1+x3-c1/c2);

fp=f(xp);

end

if f2>fp

xmin=xp;fmin=f(xp);

else

xmin=x2;fmin=f(x2);

end

end

結果：

clear all;f=x+20/x;[xmin,fmin]=main(f,2,10,0.01)

xmin =

4.4869

fmin =

8.9443

2-11

2-11

我的：

2-11

function [k ender]=steepest(f,x,e)

%梯度下降法,f為目標函數(兩變量x1和x2)，x為初始點,如[3;4]

syms x1 x2 m; %m為學習率

d=-[diff(f,x1);diff(f,x2)]; %分別求x1和x2的偏導數，即下降的方向

flag=1; %循環標志

k=0; %迭代次數

while(flag)

d_temp=subs(d,x1,x(1)); %將起始點代入，求得當次下降x1梯度值

d_temp=subs(d_temp,x2,x(2)); %將起始點代入，求得當次下降x2梯度值

nor=norm(d_temp); %范數

if(nor>=e)

x_temp=x+m*d_temp; %改變初始點x的值

f_temp=subs(f,x1,x_temp(1)); %將改變后的x1和x2代入目標函數

f_temp=subs(f_temp,x2,x_temp(2));

h=diff(f_temp,m); %對m求導，找出最佳學習率

m_temp=solve(h); %求方程，得到當次m

x=x+m_temp*d_temp; %更新起始點x

k=k+1;

else

flag=0;

end

end

ender=double(x); %終點

end

syms x1 x2;

f=x1^2+x2^2-x1*x2-10*x1-4*x2+60;

x=[0;0];

e=0.01;

[k ender]=steepest(f,x,e)

ender =

7.9961

5.9971

室友的：

2-11:

梯度函數：

function [k,ender]=tidu(f,x,e)

syms x1 x2 m;

d=-[diff(f,x1);diff(f,x2)];

flag=1;

k=0;

while(flag)

d_temp=subs(d,x1,x(1));

d_temp=subs(d_temp,x2,x(2));

nor=norm(d_temp);

if(nor>=e)

x_temp=x+m*d_temp;

f_temp=subs(f,x1,x_temp(1));

f_temp=subs(f_temp,x2,x_temp(2));

h=diff(f_temp,m);

m_temp=solve(h);

x=x+m_temp*d_temp;

k=k+1;

else

flag=0;

end

end

ender=double(x);

end

結果指令：

syms x1 x2;

f=x1^2+x2^2-x1*x2-10*x1-4*x2+60;

x=[0;0];

e=0.01;

[k ender]=tidu(f,x,e)

ender =

7.9961

5.9971

2-12

我的：

2-12

展開為二階泰勒式

syms x1 x2;

taylor(x1^4+2*x2^3-3*x1^2*x2)

ans =

3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) + 3*(x1 - 1)^2 + 6*(x2 - 1)^2 - 1

牛頓法求解：

function all=newton(f,x,e)

syms x1 x2 h;

d=-[diff(f,x1);diff(f,x2)];

h=hessian(f);

flag=1;

h1=h^-1;

while (flag)

d_temp=subs(d,x1,x(1));

d_temp=subs(d_temp,x2,x(2));

nor=norm(d_temp);

if(nor>=e)

x=x+h1*d_temp;

else

flag=0;

end

end

all=double(x);

結果指令：

clear all

>> syms x1 x2;

f=3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) + 3*(x1 - 1)^2 + 6*(x2 - 1)^2 - 1;

x=[1;1];

e=0.01;all=newton(f,x,e)

all =

1.1667

0.8333

室友的：

2-12:

展開為二階泰勒式

syms x1 x2;

taylor(x1^4+2*x2^3-3*x1^2*x2)

ans =

3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) + 3*(x1 - 1)^2 + 6*(x2 - 1)^2 - 1

牛頓函數：

function all=newton(f,x,e)

syms x1 x2 h;

d=-[diff(f,x1);diff(f,x2)];

h=hessian(f);

flag=1;

h1=h^-1;

while (flag)

d_temp=subs(d,x1,x(1));

d_temp=subs(d_temp,x2,x(2));

nor=norm(d_temp);

if(nor>=e)

x=x+h1*d_temp;

else

flag=0;

end

end

all=double(x);

結果指令：

clear all

>> syms x1 x2;

f=3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) + 3*(x1 - 1)^2 + 6*(x2 - 1)^2 - 1;

x=[1;1];

e=0.01;all=newton(f,x,e)

all =

1.1667

0.8333

2-13(1)

我的：

2-13(1)

外點懲罰函數法：

function [ x,y ] = Epfm_min( fx,gx,hx,xx0,s,c,a)

%fx是目標函數

%gx是不等式約束方程組(且g>=0)

%xx0是初始點

%hx是等式約束方程組(且h=0)

%s是精確度(s>0)

%c是放大系數(c>1)

%a是罰因子(默認為1)

syms x1 x2

xx1=xx0;

v=[x1,x2];

a1=a;

Pxk=1;%假設Px等于1，以免不必要錯誤

G=-subs(gx,v,xx1);%用于判別max{0,-g(x)}

while Pxk>s

if(G<0)

Px=a1*hx*hx;

else

Px=a1*hx*hx+a1*gx*gx;

end

Fx=fx+a1*Px;%將約束問題化為了一個無約束的問題

% 接下來解min F(x)

dFx1=diff(Fx,x1);%分別對x1，x2求偏導數

dFx2=diff(Fx,x2);

[k,b]=solve(dFx1,dFx2,'x1','x2');%求出

xx2=xx1+[k,b];

Pxk=a1*subs(Px,v,xx2);

xx1=xx2;%相當于置k=k+1

a1=c*a1;%罰因子放大

G=-subs(gx,v,xx1);%用于判別max{0,-g(x)}

end

x=xx1;

y=a1/c;

syms x1 x2;

fx=x1+x2;

gx=-x1;

hx=x1^2-x2;

s=10.^-5

c=10

xx0=[0,0]

a=1;

[x,y]=Epfm_min( fx,gx,hx,xx0,s,c,a)

>> x=[0.1;0.2];k=0.1;e=0.01;r=1;[x,minf]= Epfm_min (p,x,k,r,e)

x =

0.0015

minf =

0.0030

Ans=

0.0045

內點懲罰函數法

function [x,minf]=minNF(f,x0,g,u,v,var,eps)

format long;

if nargin==6

eps=1.0e-4;

end

k=0;

FE=0;

for i=1:length(g)

FE=FE+1/g(i);

end

x1=transpose(x0);

x2=inf;

while 1

FF=u*FE;

SumF=f+FF;

[x2,minf]=minNT(SumF,transpose(x1),var);

Bx=Funval(FE,var,x2);

if u*Bx

if norm(x2-x1)<=eps

x=x2;

break;

else

u=v*u;

x1=x2;

end

else

if norm(x2-x1)<=eps

x=x2;

break;

else

u=v*u;

x1=x2;

end

end

end

minf=Funval(f,var,x);

format short;

syms x1 x2 r1;

>> p=taylor(x1+x2-r1*(1/(x1^2-x2)-1/x1),[x1 x2],[0.001 0.002],'Order',3);

>> x=[0.1;0.2];k=0.1;e=0.01;r=1;[x,minf]= minNF(p,x,k,r,e)

x =

0.0015

minf =

0.0030

Ans=

0.0045

室友的：

2-13：

內點：

懲罰函數：

function anll=neicheng(p,x,k,r,e)

syms x1 x2 r1;

flag1=1;

while (flag1)

pd=subs(p,r1,r);

xold=x;

flag2=1;

while (flag2)

dp=-[diff(pd,x1);diff(pd,x2)];

h=hessian(pd,[x1,x2]);

h1=h^-1;

dp_temp=subs(dp,x1,x(1));

dp_temp=subs(dp_temp,x2,x(2));

nor=norm(dp_temp);

if(nor>=e)

x=x+h1*dp_temp;

else

flag2=0;

end

end

x_temp=x;

nor2=norm(x_temp-xold);

if double(nor2)>=e

r=k*r;

else

flag1=0;

end

end

anll=double(x);

結果：

clear all; syms x1 x2 r1;

>> p=taylor(x1+x2-r1*(1/(x1^2-x2)-1/x1),[x1 x2],[0.001 0.002],'Order',3);

>> x=[0.1;0.2];k=0.1;e=0.01;r=1;anll=neicheng(p,x,k,r,e)

anll =

0.0015

0.0030

Ans=

0.0045

外點：

懲罰函數：

function annn=waicheng(p,x,k,r,e)

syms x1 x2 r1;

flag1=1;

while (flag1)

pd=subs(p,r1,r);

xold=x;

flag2=1;

while (flag2)

dp=-[diff(pd,x1);diff(pd,x2)];

h=hessian(pd,[x1,x2]);

h1=h^-1;

dp_temp=subs(dp,x1,x(1));

dp_temp=subs(dp_temp,x2,x(2));

nor=norm(dp_temp);

if(nor>=e)

x=x+h1*dp_temp;

else

flag2=0;

end

end

x_temp=x;

nor2=norm(x_temp-xold);

if double(nor2)>=e

r=k*r;

else

flag1=0;

end

end

annn=double(x);

結果：

clear all; syms x1 x2 r1;

p=taylor(x1+x2,[x1 x2],[0.001 0.002],'Order',3);

x=[0.1;0.2];k=0.1;e=0.01;r=1;annn=waicheng(p,x,k,r,e)

annn =

0.0015

0.0030

Ans=

0.0045

2-14

我的(貌似這題抄的他的)：

2-14

function [x,minf] = minMixFun(f,g,h,x0,r0,c,var,eps)

gx0 = Funval(g,var,x0);

if gx0 >= 0;

else

disp('初始點必須滿足不等式約束！');

x = NaN;

minf = NaN;

return;

end

if r0 <= 0

disp('初始障礙因子必須大于0！');

x = NaN;

minf = NaN;

return;

end

if c >= 1 || c < 0

disp('縮小系數必須大于0且小于1！');

x = NaN;

minf = NaN;

return;

end

if nargin == 7

eps = 1.0e-6;

end

FE = 0;

for i=1:length(g)

FE = FE + 1/g(i);

end

FH = transpose(h)*h;

x1 = transpose(x0);

x2 = inf;

while 1

FF = r0*FE + FH/sqrt(r0);

SumF = f + FF ;

[x2,minf] = minNT(SumF,transpose(x1),var);

if norm(x2 - x1)<=eps

x = x2;

break;

else

r0 = c*r0;

x1 = x2;

end

end

minf = Funval(f,var,x);

Funval.m

function fv = Funval(f,varvec,varval)

var = findsym(f);

varc = findsym(varvec);

s1 = length(var);

s2 = length(varc);

m =floor((s1-1)/3+1);

varv = zeros(1,m);

if s1 ~= s2

for i=0: ((s1-1)/3)

k = findstr(varc,var(3*i+1));

index = (k-1)/3;

varv(i+1) = varval(index+1);

end

fv = subs(f,var,varv);

else

fv = subs(f,varvec,varval);

end

Syms x1 x2;

f=x1^2-x2^2-3*x2;

g=1-x1;

h=x2-2;

[x,minf]=minMixFun(f,g,h,[2,2],2,0.5,[x1 x2 ],0.001)

x =

1.0015

minf=

2.0002

室友的：

2-14：

混合：

懲罰函數：

function [x,minf] = MixPunish(f,g,h,x0,r0,c,var,eps)

gx0 = Funval(g,var,x0);

if gx0 >= 0;

else

disp('初始點必須滿足不等式約束！');

x = NaN;

minf = NaN;

return;

end

if r0 <= 0

disp('初始障礙因子必須大于0！');

x = NaN;

minf = NaN;

return;

end

if c >= 1 || c < 0

disp('縮小系數必須大于0且小于1！');

x = NaN;

minf = NaN;

return;

end

if nargin == 7

eps = 1.0e-6;

end

FE = 0;

for i=1:length(g)

FE = FE + 1/g(i);

end

FH = transpose(h)*h;

x1 = transpose(x0);

x2 = inf;

while 1

FF = r0*FE + FH/sqrt(r0);

SumF = f + FF ;

[x2,minf] = minNT(SumF,transpose(x1),var);

if norm(x2 - x1)<=eps

x = x2;

break;

else

r0 = c*r0;

x1 = x2;

end

end

minf = Funval(f,var,x);

Funval.m

function fv = Funval(f,varvec,varval)

var = findsym(f);

varc = findsym(varvec);

s1 = length(var);

s2 = length(varc);

m =floor((s1-1)/3+1);

varv = zeros(1,m);

if s1 ~= s2

for i=0: ((s1-1)/3)

k = findstr(varc,var(3*i+1));

index = (k-1)/3;

varv(i+1) = varval(index+1);

end

fv = subs(f,var,varv);

else

fv = subs(f,varvec,varval);

end

Syms x1 x2;

f=x1^2-x2^2-3*x2;

g=1-x1;

h=x2-2;

[x,minf]=MixPunish(f,g,h,[2,2],2,0.5,[x1 x2 ],0.001)

x =

1.0015

minf=

2.0002

結束語

無聊到這地步想必也是沒誰了。不過，剛考完，我總不能一直玩手機啊。前幾天重新看《盤龍》讓我在手機上剛了一星期，不能再這么毫無節制的玩耍了，但是又不想學習，所以只好寫簡書了~~

不過網上畢竟這方面的資源不是很多，我就算是為后來人做點好事吧，讓你們好找一點~~

個人宣言

知識傳遞力量，技術無國界，文化改變生活！