注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

信息 灵感 创新

I? =Information,Inspiration,Innovation

 
 
 

日志

 
 
关于我

we are 5. Mathematics, Computation, Programming, Engineering, and Making fun of life.

网易考拉推荐

【转载】matlab的一些函数模板与巧妙的方法  

2013-01-02 01:57:03|  分类: M&M |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
最近在看《matlab高效编程25个案例分析》一书,里面提到的一些算法确实很好,于是就做了如下的摘录:
%% application of anonymous function
%1.build a anonymous function
f=@(x,y)2*x+y.^3;               %constuct the handle of function
f(1,2);
f=@(a,b)@(x)a*x+b;
f23=f(2,3);
f23info=functions(f23);

%2.equation solve
f=@(a)@(x)exp(x)+x^a+x^(sqrt(x))-100;
format long;
aa=0:.01:2;
plot(aa,arrayfun(@(a)fzero(f(a),4),aa),'* -');   %apply the use of arrayfun
xlabel('$a$','interpreter','latex','fontsize',15);
ylabel('$x$','interpreter','latex','fontsize',15);
title('$\mathrm{e}^{x}+x^{\sqrt{x}}+x^a-100$','interpreter','latex','fontsize',15);  %label with latex format

%3.describe implicit function
y=@(x)fzero(@(y)(exp(y)+x^y)^(1/y)-x^2*y,1);
y=@(x)arrayfun(@(xx)fzero(@(y)(exp(y)+xx^y)^(1/y)-xx^2*y,1),x);  %vectorized function

%4.differential problems
syms x
f=(x+tan(x))^(sin(x));
c=diff(f,3);
f3=eval(['@(x)' vectorize(c)]);
x=linspace(0,1,100);
plot(x,f3(x),'LineWidth',2)
title('y=(x+tan(x))\^(sin(x))三阶导数的图像');
xlabel('x');ylabel('y');

%% application of integral
%1.comparasion of three kinds of methods for triple integral
tic;quadl(@(x)arrayfun(@(xx)quad2d(@(y,z)xx.*y.*z,xx,2*xx,@(y)xx*y,@(y)...
    2*xx*y),x),1,2),toc;
tic;quad2d(@(x,y)arrayfun(@(xx,yy)quadl(@(z)...
    xx.*yy.*z,xx*yy,2*xx*yy),x,y),1,2,@(x)x,@(x)2*x),toc;
tic;quadl(@(x)arrayfun(@(xx)quadl(@(y)arrayfun(@(yy)quadl(@(z)...
    xx.*yy.*z,xx*yy,2*xx*yy),y),xx,2*xx),x),1,2),toc;
tic;triplequad(@(x,y,z)x.*y.*z.*(z<=2*x.*y&z>=x.*y&y<=2*x&y>=x),...
    1,2,1,4,1,16),toc;
%notes:It is obvious that the combination of basic functions would run mush
%more faster

%% qpplication of optimazation
%1.randwalk method and its illustration
f2=@(x,y)-sin(sqrt((x-50).^2+(y-50).^2)+exp(1))./(sqrt((x-50).^2+(y-50).^2)+exp(1))-1;
[mx,minf]=randwalknew(f2,[0 0],10,0.00001,1000,10);
% an example of implicit function
z=@(x,y)fzero(@(z)z-sin((z*x-0.5)^2+x*2*y^2-z/10)*exp(-((x-0.5-exp(-y+z))^2+...
    y^2-z/5+3)),rand);
tic,[mx,minf]=randwalknew(z,[unifrnd(-1,7),unifrnd(-2,2)],5,0.000001,100,5),toc;

%2.find the solve of nonlinear equation group
p=fsolveDemo1;
disp(num2str(p));

%3.find the solution of equation with integral
ker=@(x,s)x*s-1;
f=@(x)polyval([-1/4,-2/3,-1/6,1,2,1],x);%自由项函数
lambda=1;
FY=@(Y)Y^2;%y(x)的非线性泛函
a=0;
b=1;
n=10;%10等分[a,b]区间
Y0=ones(1,n+1);%求解Y的初始值
[X,Y,fval]=NonlinearVolterra2(ker,f,lambda,FY,a,b,n,Y0);

%4.find the solution of mixture binary optimazation
 C=[5 -6];
 D=[7 3];
 A=[-2 3;2 -1;1 1];
 B=[1 0;1 1;0 1];
 b=[0;11;7];
 [OptX,OptY,OptValue]=BendersDecomposition(C,D,A,B,b);
关于以上代码中提到的自定义函数如下:

首先是BendersDecomposition函数:
function [OptX,OptY,OptValue]=BendersDecomposition(C,D,A,B,b)
%用于解决如下的问题:
%       min    C*x+D*y
%       s.t.     A*x+B*y>=b;x,0-1变量,y>=0

%迭代停止的判断阀值:LB/UB大于epsilon时停止迭代。
epsilon=0.99;
[nRowA,nColA]=size(A);
[nRowB,nColB]=size(B);
%step1:初值化
x0=zeros(nColA,1);%初始值
LB=-1e10;
UB=inf;
p=0;
q=0;
U=zeros(100,nRowA);
V=zeros(100,nRowA);
options=optimset('LargeScale','off','Simplex','on');
OptionsBint=optimset('MaxRLPIter',100000,'NodeSearchStrategy','bn',...
    'MaxTime',50000);
while LB/UB<epsilon
    %step2:求解子问题
    [UorV,fval,exitflag]=linprog((A*x0-b),B',D',[],[],zeros(nRowB,1),...
        inf(nRowB,1),[],options);
    if exitflag==1
        p=p+1;
        U(p,:)=UorV';%极值点
        UB=C*x0+U(p,:)*(b-A*x0);
    elseif exitflag==-3 %无界
        q=q+1;
        V(q,:)=UorV'/1e16;%极方向
        V(q,V(q,:)<0.0001)=0;
    else
        msgbox('无可行解!');
        return;
    end
    %step3:求解主问题
    if isinf(UB)
        k=17;
    else
        k=ceil(log2(UB+1))-1;%根据2^(k+1)-1>UB决定最小的k值
    end
    CoefZ=2.^(0:k); %z0到zk的系数向量
    f=[CoefZ';zeros(nColA,1)];%主问题目标函数的系数向量
    CoefMatZX=[repmat(-CoefZ,p,1),-U(1:p,:)*A;zeros(q,k+1),-V(1:q,:)*A];
    bZX=[-(repmat(C*x0,p,1)+U(1:p,:)*b);-V(1:q,:)*b];
    [OptZX,minZ,ExitflagBint]=bintprog(f,CoefMatZX,bZX,[],[],[],...
        OptionsBint);
        if ExitflagBint==1
            LB=minZ;
            x0=OptZX(k+2:end);
        else
            break;
        end
        LB,UB % #OK<NOPRT>
end
OptX=x0;
options=optimset('LargeScale','off','Simplex','on');
[OptY,OptValue]=linprog(D',-B,A*OptX-b,[],[],zeros(nColB,1),inf(nColB,1),...
    [],options);
OptValue=C*OptX+D*OptY;


接下来是fsolveDemo1:
function p=fsolveDemo1
p0=[51.61;-1;0;1];
options=optimset('MaxFunEvals',20000,'MaxIter',2000);
[p,fval]=fsolve(@f,p0,options);
    function F=f(p) %待求方程的解
        F=[p(1)+p(2)*(1-exp(-(p(3)*(0)^p(4))))-51.61;
            p(1)+p(2)*(1-exp(-p(3)*(9.78)^p(4)))-51.91;
            p(1)+p(2)*(1-exp(-(p(3)*(30.68)^p(4))))-53.27;
            p(1)+p(2)*(1-exp(-(p(3)*(59.7)^p(4))))-59.68;];
    end
end

然后是NonlinearVolterra2函数:
function [X,Y,fval]=NonlinearVolterra2(ker,f,lambda,FY,a,b,n,Y0)
%ker:核函数的函数句柄
%f:自由项f(x)的函数句柄
%lambda:参数
%a,b:积分u区间
%n:积分区间被等分的个数
%Y0:求解Y的时候初值向量
if length(Y0)~=n+1
    error(['Y0的长度必须要等于n+1:',num2str(n+1)]);
end
Y0=Y0(:);
eval(['ker=',vectorize(ker),';']);%ker为统一的点乘形式
eval(['f=',vectorize(f),';']);%f同上
eval(['FY=',vectorize(FY),';']);%FY同上
if nargin<=5%如果没有指定n的值
    n=5;
end
A=cell(n,1);%存储cotes系数
for k=1:n
    A{k}=NewtonCotes(k);
end
K=zeros(n+1);
h=(b-a)/n;
xk=(a:h:b)';
for ix=2:n+1
    for jx=1:ix
        K(ix,jx)=(ix-1)*h*A{ix-1}(jx)*ker(xk(ix),xk(jx));%计算k矩阵
    end
end
K=lambda*K;
F=f(xk);
if isscalar(F)%针对某些自由函数为常数的情形
    F=F*ones(size(xk));
end
X=xk';
options=optimset('TolFun',1e-16);%非线性方程求解设置
%利用fsolve函数求解离散原方程后的到非线性方程组
[Y,fval]=fsolve(@(Y)K*FY(Y)+F-Y,Y0,options);
Y=Y';
fval=fval';
end
%===================================
%求cotes系数的函数(普通子函数形式)
%===================================
function A=NewtonCotes(n)
%n等分cotes区间,返回cotes系数
syms t;
A=zeros(n+1,1);
for k=0:n
    %利用符号计算得到Cotes系数中被积函数的表达式
    symfun=prod(t-[0:k-1,k+1:n]);
    %将得到的被积函数表达式转成“点乘”的形式并生成相应的函数句柄
    eval(['fun=@(t)',vectorize(symfun),';']);
    A(k+1)=quadgk(fun,0,n)*(-1)^(n-k)/gamma(k+1)/gamma(n-k+1)/n;
end
end





最后是随机行走函数randwalknew
function [mx,minf]=randwalknew(f,x,lamda,epsilon,N,n)
%f is the handle of the function
%x is the initial value.lamda acts as steps
%eplision is the Threshold values  of the precision
%N is the num of the loops,n is thenum of steps
F=zeros(n,1);
X=zeros(n,2);
f1=f(x(1),x(2));
while(lamda>=epsilon)
    k=1;
    while(k<=N)
        u=unifrnd(-1,1,n,2);
        for ii=1:n
            X(ii,:)=x+lamda*(u(ii,:)/norm(u(ii,:)));
            F(ii)=f(X(ii,1),X(ii,2));
        end
        [f11,kk]=min(F);
        if f11<f1
            f1=f11;
            x=X(kk,:);
            k=1;
        else
            k=k+1;
        end
    end
    lamda=lamda/2;
end
mx=X(kk,:);
minf=f1;
end
%% the end of the randwalk function


  评论这张
 
阅读(854)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016