Advanced Mathematics and Mechanics Applications Using MATLAB phần 10 pdf

148: num2str(dbest),’, CPU TIME = ’,...

149: num2str(t),’ SECS’])

150: rotate3d on, shg, disp(’ ’) 151: disp(’Rotate the figure or press’)

152: disp(’return to continue’) 153: dumy=input(’ ’,’s’); close

154:

155: end 156:

157: %===========================================

158:

159: function r=cylpoint(w1,w2,r0,m,rdat,zdat)

160: % r=cylpoint(w1,w2,v,r0,m,rdat,zdat)

161: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 162: % This function computes the position of a

163: % point on the surface of a circular cylinder 164: % arbitrarily positioned in space. The argument

165: % list parameters have the following form,

166: % where rad means cylinder radius, and len 167: % means cylinder length.

168: % b=2*rad+len;

169: % zdat=[[0,0]; [rad/b, 0]; 170: % [(rad+len)/b, len];[ 1, len]];

171: % rdat=zdat; rdat(2,2)=rad; 172: % rdat(3,2)=rad; rdat(4,2)=0;

173:

174: u=2*pi*sin(w1)^2; v=sin(w2)^2; 175: z=interp1(zdat(:,1),zdat(:,2),v);

176: rho=interp1(rdat(:,1),rdat(:,2),v);

177: x=rho*cos(u); y=rho*sin(u); 178: r=r0(:)+m*[x;y;z];

179:

180: %===========================================

181:

182: function dsqr=dcyl2cyl(... 183: w,r0,m,rdat,zdat,R0,M,Rdat,Zdat)

184: %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

185: % dsqr=dcyl2cyl(w,r0,m,rdat,zdat,R0,M,Rdat,Zdat) 186: % This function computes the square of the

187: % distance between generic points on the 188: % surfaces of two circular cylinders in three

189: % dimensions.

190: % 191: % User m functions called: cylpoint

192:

© 2003 by CRC Press LLC

193: global fcount

194: fcount=fcount+1;

195: r=cylpoint(w(1),w(2),r0,m,rdat,zdat); 196: R=cylpoint(w(3),w(4),R0,M,Rdat,Zdat);

197: dsqr=norm(r-R)^2; 198:

199: %===========================================

200:

201: function cylfigs

202: % cylfigs

203: % ~~~~~~~ 204: % This function plots the geometries

205: % pertaining to four data cases used

206: % to test closest proximity problems 207: % involving two circular cylinders

208: % 209: % User m functions called: plot2cyls

210:

211: w=rads; p=1:2; q=3:4; s=5:6; t=7:8; 212:

213: rad=1; len=3; r0=[4,0,0]; v=[0,0,1];

214: Rad=1; Len=3; R0=[0,4,0]; V=[0,0,1]; 215: d=.4; subplot(2,2,1)

216: [x,y,z,X,Y,Z]=plot2cyls(rad,len,r0,v,Rad,Len,... 217: R0,V,d,’CASE 1’); hold on

218: plot3(w(p,1),w(p,2),w(p,3),’linewidth’,2’)

219: hold off 220:

221: rad=1; len=3; r0=[4,0,0]; v=[3,0,4];

222: Rad=1; Len=3; R0=[0,4,0]; V=[0,3,4]; 223: d=.4; subplot(2,2,2);

224: [x,y,z,X,Y,Z]=plot2cyls(rad,len,r0,v,Rad,Len,... 225: R0,V,d,’CASE 2’); hold on

226: plot3(w(q,1),w(q,2),w(q,3),’linewidth’,2’)

227: hold off 228:

229: rad=1; len=5; r0=[4,0,0]; v=[-4,0,3];

230: Rad=1; Len=5; R0=[0,4,0]; V=[0,0,1]; 231: d=.4; subplot(2,2,3)

232: [x,y,z,X,Y,Z]=plot2cyls(rad,len,r0,v,Rad,Len,... 233: R0,V,d,’CASE 3’); hold on

234: plot3(w(s,1),w(s,2),w(s,3),’linewidth’,2’)

235: hold off 236:

237: rad=1; len=4*sqrt(2); r0=[4,0,0]; v=[-1,1,0];

© 2003 by CRC Press LLC

238: Rad=1; Len=3; R0=[0,0,-2]; V=[0,0,-1];

239: d=.4; subplot(2,2,4);

240: [x,y,z,X,Y,Z]=plot2cyls(rad,len,r0,v,Rad,Len,... 241: R0,V,d,’CASE 4’); hold on

242: plot3(w(t,1),w(t,2),w(t,3),’linewidth’,2’) 243: hold off, subplot

244: % print -deps cylclose

245:

246: %===========================================

247:

248: function [x,y,z,X,Y,Z]=plot2cyls(... 249: rad,len,r0,vc,Rad,Len,R0,Vc,d,titl)

250: % [x,y,z,X,Y,Z]=plot2cyls(rad,len,r0,vc,Rad,...

251: % Len,R0,Vc,d,titl) 252: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

253: % This function generates point grids on the 254: % surfaces of two circular cylinders and plots

255: % both cylinders together

256: % 257: % User m functions called: cornrpts surfmany

258: % cylpts

259: if nargin==0 260: titl=’TWO CYLINDERS’;

261: rad=1; len=3; r0=[4,0,0]; vc=[3,0,4]; 262: Rad=1; Len=3; R0=[0,4,0]; Vc=[0,3,4]; d=.2;

263: end

264: if isempty(titl), titl=’ ’; end 265: u=2*rad+len; v=2*pi*rad;

266: nu=ceil(u/d); nv=ceil(v/d);

267: u=cornrpts([0,rad,rad+len,u],nu)/u; 268: v=linspace(0,1,nv);

269: [x,y,z]=cylpts(u,v,rad,len,r0,vc); 270: U=2*Rad+Len; V=2*pi*Rad;

271: Nu=ceil(U/d); Nv=ceil(V/d);

272: U=cornrpts([0,Rad,Rad+Len,U],Nu)/U; 273: V=linspace(0,1,Nv);

274: [X,Y,Z]=cylpts(U,V,Rad,Len,R0,Vc);

275: surfmany(x,y,z,X,Y,Z), title(titl) 276: colormap([1 1 0]), shg

277:

278: %===========================================

279:

280: function [x,y,z]=cylpts(... 281: axial,circum,rad,len,r0,vectax)

282: % [x,y,z]=cylpts(axial,circum,rad,len,r0,vectax)

© 2003 by CRC Press LLC

283: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

284: % This function computes a grid of points on the

285: % surface of a circular cylinder 286: %

287: % User m functions called: ortbasis 288:

289: U=2*rad+len; u=U*axial(:); n=length(u);

290: v=2*pi*circum(:)’; m=length(v); 291: ud=[0,rad,rad+len,U];

292: r=interp1(ud,[0,rad,rad,0],u);

293: z=interp1(ud,[0,0,len,len],u); 294: x=r*cos(v); y=r*sin(v); z=repmat(z,1,m);

295: % w=basis(vectax)*[x(:),y(:),z(:)]’;

296: w=ortbasis(vectax)*[x(:),y(:),z(:)]’; 297:

298: x=r0(1)+reshape(w(1,:),n,m); 299: y=r0(2)+reshape(w(2,:),n,m);

300: z=r0(3)+reshape(w(3,:),n,m);

301:

302: %===========================================

303:

304: function v=cornrpts(u,N) 305: % v=cornrpts(u,N)

306: % ~~~~~~~~~~~~~~ 307: % This function generates approximately N

308: % points between min(u) and max(u) including

309: % all points in u plus additional points evenly 310: % spaced in each successive interval.

311: % u - vector of points

312: % N - approximate number of output points 313: % between min(u(:)) and max(u(:))

314: % v - vector of points in increasing order 315:

316: u=sort(u(:))’; np=length(u);

317: d=u(np)-u(1); v=u(1); 318: for j=1:np-1

319: dj=u(j+1)-u(j); nj=max(1,fix(N*dj/d));

320: v=[v,[u(j)+dj/nj*(1:nj)]]; 321: end

322:

323: %===========================================

324:

325: function mat=ortbasis(v) 326: % mat=ortbasis(v)

327: % ~~~~~~~~~~~~~~

© 2003 by CRC Press LLC

328: % This function generates a rotation matrix

329: % having v(:)/norm(v) as the third column

330:

331: v=v(:)/norm(v); mat=[null(v’),v];

332: if det(mat)<0, mat(:,1)=-mat(:,1); end 333:

334: %===========================================

335:

336: function [xmin,fmin,m,ntype]=nelmed(...

337: F,x0,dx,epsx,epsf,M,ifpr,varargin)

338: % [xmin,fmin,m,ntype]=nelmed(... 339: % F,x0,dx,epsx,epsf,M,ifpr,varargin)

340: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

341: % This function performs multidimensional 342: % unconstrained function minimization using the

343: % direct search procedure developed by 344: % J. A. Nelder and R. Mead. The method is

345: % described in various books such as:

346: % ’Nonlinear Optimization’, by M. Avriel 347: %

348: % F - objective function of the form

349: % F(x,p1,p2,...) where x is vector 350: % in n space and p1,p2,... are any

351: % auxiliary parameters needed to 352: % define F

353: % x0 - starting vector to initiate

354: % the search 355: % dx - initial polyhedron side length

356: % epsx - convergence tolerance on x

357: % epsf - convergence tolerance on 358: % function values

359: % M - function evaluation limit to 360: % terminate search

361: % ifpr - when this parameter equals one,

362: % different stages in the search 363: % are printed

364: % varargin - variable length list of parameters

365: % which can be passed to function F 366: % xmin - coordinates of the smallest

367: % function value 368: % fmin - smallest function value found

369: % m - total number of function

370: % evaluations made 371: % ntype - a vector containing

372: % [ninit,nrefl,nexpn,ncontr,nshrnk]

© 2003 by CRC Press LLC

373: % which tells the number of reflect￾374: % ions, expansions, contractions,and

375: % shrinkages performed 376: %

377: % User m functions called: objective function 378: % named in the argument list

379:

380: if isempty(ifpr), ifpr=0; end 381: if isempty(M), M=500; end;

382: if isempty(epsf), epsf=1e-5; end

383: if isempty(epsx), epsx=1e-5; end 384:

385: % Initialize the simplex array

386: x0=x0(:); n=length(x0); N=n+1; f=zeros(1,N); 387: x=repmat(x0,1,N)+[zeros(n,1),dx*eye(n,n)];

388: for k=1:N 389: f(k)=feval(F,x(:,k),varargin{:});

390: end

391:

392: ninit=N; nrefl=0; nexpn=0; ncontr=0;

393: nshrnk=0; m=N;

394:

395: Erx=realmax; Erf=realmax;

396: alpha=1.0; % Reflection coefficient 397: beta= 0.5; % Contraction coefficient

398: gamma=2.0; % Expansion coefficient

399:

400: % Top of the minimization loop

401:

402: while Erx>epsx | Erf>epsf 403:

404: [f,k]=sort(f); x=x(:,k); 405:

406: % Exit if maximum allowable number of

407: % function values is exceeded 408: if m>M, xmin=x(:,1); fmin=f(1); return; end

409:

410: % Generate the reflected point and 411: % function value

412: c=sum(x(:,1:n),2)/n; xr=c+alpha*(c-x(:,N)); 413: fr=feval(F,xr,varargin{:}); m=m+1;

414: nrefl=nrefl+1;

415: if ifpr==1, fprintf(’ :RFL \n’); end 416:

417: if fr<f(1)

© 2003 by CRC Press LLC