# to unbundle, sh this file (in an empty directory) echo README 1>&2 sed >README <<'//GO.SYSIN DD README' 's/^-//' - - COCA User's Guide - - Description of a collection of Matlab function files - for the computation of - linear Chebyshev approximations in the complex plane - - - Bernd Fischer and Jan Modersitzki - - University of Hamburg - Institute of Applied Mathematics - Bundesstrasse 55 - D - 2000 Hamburg 13, F.R.G. - - Fax : +40-4123-5117 - e-mail: na.fischer@na-net.ornl.gov - na.modersitzki@na-net.ornl.gov - - - -1. Introduction - -COCA (COmplex linear Chebyshev Approxiamtion) is a collection -of Matlab function files for computing the best Chebyshev approximation -to a function F on a region \Omega (with boundary \partial \Omega) -with respect to the n-dimensional linear function space V -spanned by the functions p_1 ,p_2 ,...,p_n : - - max | F(z) - sum_{j=1}^n a_j p_j(z)| = Min! - z \in \Omega - -It is assumed that the maximum principle holds, i.e., the Chebyshev -norm over \Omega can be replaced by the one over the boundary \partial -\Omega. - -The mathematical background of the algorithms of COCA will be -described in a forthcoming paper. - - -2. How to obtain COCA - -It is planned to submit the program to the Matlab part of the "Netlib" -facility. Until then, please contact one of the authors. - -COCA only uses standard Matlab files from the general Matlab -toolbox, so it should run in any Matlab environment. - - -3. Use of COCA - -COCA (including 17 examples and 17 predefined functions, bases, -and boundaries) consists of 56 (sorry about this) m-files, one -of which is most important for the user - - COCA: drives the following two main ingredients of the package - - REMEZ: basically solves the problem - Q_NEWTON: refines results of REMEZ - -Ouit interesting is in addition the file - - SET_PARA: sets various parameters - -cause this is the place to learn about the default settings of certain -parameters which affect the performance of the iterative algorithms. -In addition, the level of output (including plots) is controlled -by these parameters. - -Here are all routines of COCA: - -BOUNDS DEF_VAR NEWTON_F SET_COLM STOPITER -CLUSTER ERRFUN PICTURE SET_EXTR STOPMENU -COCA (*) IMPROVE PLOT_EXT SET_PARA (*) -COMPNORM INITIAL Q_NEWTON (*) SIGN_SPE -CONVSTR MENUSTR REMEZ (*) SIMPLEX - - -In order to compute the best approximation to a function F, on a -boundary \partial \Omega with respect to an -approximation space V, the user has to supply three function files. - -1) A function file which defines the function F to be approximated. - -2) A function file which computes the boundary \partial \Omega in - terms of a parametrization \gamma over the interval [0,1], i.e., - \gamma([0,1]) = \partial \Omega. - -3) A function file which evaluates the basis functions in V for any - given value. - -It is assumed that all three subroutine can "live" with vectors -as input, e.g., the monom z^n should be coded as z.^n. - -It is worthwhile to mention that already quite a few m-files for the -definition of various functions, bases, and boundaries are included in -the package: - -FCOS MONOM CIRCLE -FEXP MONOM_C CIRCLE_C -FMONOM ELLIPSE -FONE ELLIPS_C -FPOLE_C INTERVAL -FSIN L_SHAPE - RECANGLE - SECTOR - TWO_CIRC - -If the user has in mind to apply as well Newtons method, it is -necessary to compute in addition the corresponding derivatives. - -Finally the user has to call the function COCA. - -This is all the user has to do. However, before running the first own -problem, it is highly recommended "to play" with the examples - -EXAM_IJ.M, I=1,2,3,4,5 and J=1,2,3,(4) - -which come with the package. They are designed to show most of the -capabilities of COCA and should be easy to understand (at least -from our point of view). - - -Comments in any form are greatly appreciated. - - -Final note: This user's guide is a preliminary version of the real one, -which will be (hopefully) available in the near future. - - - -Hamburg, 17.07.92 - -Bernd Fischer, Jan Modersizki - //GO.SYSIN DD README echo bounds.m 1>&2 sed >bounds.m <<'//GO.SYSIN DD bounds.m' 's/^-//' - function[ x, error_norm, relative_error, t_max, alpha_max ] = ... - bounds( param, FUNstr, BASISstr, BOUNDARYstr, x, A, c ); -%function[ x, error_norm, relative_error, t_max, alpha_max ] = ... -%bounds( param, FUNstr, BASISstr, BOUNDARYstr, x, A, c ); -% -% BOUNDS.M 16.07.92 -% -% Computes the coefficients of the actual approximation and a lower -% and upper bound for the minimal deviation -% (for details see: COMPNORM.M). -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means, manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - R_dim = param(2); - iter = param(9); - output = param(17); - m = R_dim + 1; - - if output > 0 - disp(['#Iter ', int2str(iter)]); - end; - - if output > 1 - disp([ 'BOUNDS']) - end % if - - x( 3*m+1:4*m ) = A'\c; - - [ error_norm, relative_error, t_max, alpha_max ] = ... - compnorm( param, FUNstr, BASISstr, BOUNDARYstr, x ); -return; //GO.SYSIN DD bounds.m echo circle.m 1>&2 sed >circle.m <<'//GO.SYSIN DD circle.m' 's/^-//' - function [gamma, dgamma] = circle(t); -%function [gamma, dgamma] = circle(t); -% -% CIRCLE.M 16.07.92 -% -% Unit circle. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - t = 2*pi*t; - gamma = cos(t) + i*sin(t); - dgamma = (-sin(t) + i*cos(t))*2*pi; -return; //GO.SYSIN DD circle.m echo circle_c.m 1>&2 sed >circle_c.m <<'//GO.SYSIN DD circle_c.m' 's/^-//' - function [gamma, dgamma] = circle_c(t, Rec, Imc, r); -%function [gamma, dgamma] = circle_c(t); -% -% CIRCLE_C.M 16.07.92 -% -% Circle with centre c = Rec + i*Imc and radius r. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - c = Rec + i*Imc; - t = 2*pi*t; - gamma = r * (cos(t) + i*sin(t)) + c; - dgamma = r * ( -sin(t) + i*cos(t))*2*pi; -return; //GO.SYSIN DD circle_c.m echo cluster.m 1>&2 sed >cluster.m <<'//GO.SYSIN DD cluster.m' 's/^-//' - function [ x ] = cluster( param, BOUNDARYstr, x ); -%function [ x ] = cluster( param, BOUNDARYstr, x ); -% -% CLUSTER.M 16.07.92 -% -% This function enables a picture supported interactive clustering -% of extremal points. -% -% Two (or more) 'close' extremal points (t_k,alpha_k), k=1,2 -% are substituted by their weighted centre: -% -% (t,alpha) = (\sum_{k=1,2} r_k*(t_k,alpha_k) ) / \sum_{k=1,2} r_k. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% extract various variables from x - [length_t, t, alpha, r, lower_bound, lambda] = ... - def_var( param, x ); - - number = 1; % extremal points are numbered in PLOT_EXT.M - [t,I] = sort( t ); - alpha = alpha(I); - r = r(I); - - clc; - l = 1; - while l ~= 0 - - plot_ext( param, BOUNDARYstr, x, number) - disp( [ '[1 2 7] cluster points',... - ' with index number 1,2,7'] ); - disp( [ '[] stop clustering' ]); - disp( [ '[0] plot again ' ]); - index = input('select extremal points: [i1 i2 ...] ' ); - l = length(index); - - if l >= 2 - t(index(1)) = sum(t(index).*r(index)) / sum(r(index)); - t(index(2:l)) = []; - alpha(index(1)) = sum(alpha(index).*r(index)) / sum(r(index)); - alpha(index(2:l)) = []; - r(index(1)) = sum(r(index)); - r(index(2:l)) = []; - - x = [t; alpha; r; lower_bound; lambda]; - end; % if - - end; % while - -return; //GO.SYSIN DD cluster.m echo coca.m 1>&2 sed >coca.m <<'//GO.SYSIN DD coca.m' 's/^-//' - function [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); -%function [param, t, alpha, r, lower_bound, lambda, error_norm ] =... -%coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... -% BOUNDARYstr, BOUNDARYarg); -% -% COCA.M 16.07.92 -% -% This is the primary function of the COCA-package. It serves as -% a driver for different tasks: -% -% 1: REMEZ.M -% 2: Q_NEWTON.M -% 3: CLUSTER.M -% 4: PLOT_EXT.M -% 5: PICTURE.M -% -% For further informations see comments in the respective M-files. -% INPUT and OUTPUT parameters are explained in REMEZ.M and SET_PARA.M. -% -% Copyright (C) 1992, Bernd Fischer, Jan Modersitzki. -% All rights reserved. - -% -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -% convert-strings and merge with arguments - [ FUNstr, BASISstr, BOUNDARYstr ]=... - convstr( FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg, param(17) ); - -% reset variable - x = zeros(4*(param(2) + 1),1); - -% event loop: -% param(16) controls the program, -% param(15) controls the interrups - - while param(16) ~= 0 - - if param(16) == 1 - - disp( ['CALL REMEZ'] ); - - [ param, x, error_norm, relative_error ] = ... - remez( param, FUNstr, BASISstr, BOUNDARYstr, x ); - - elseif param(16) == 2 - - disp( ['CALL NEWTON'] ); - - [ param, x, error_norm, relative_error ] = ... - q_newton( param, FUNstr, BASISstr, BOUNDARYstr, x ); - - elseif param(16) == 3 - - [ x ] = cluster( param, BOUNDARYstr, x ); - param(15) = 8; % this implies a Newton call - - elseif param(16) == 4 - - number = 0; % extremal points are not numbered - plot_ext( param, BOUNDARYstr, x, number ); - - elseif param(16) == 5 - - m = param(2) + 1; - t = x( 1: m); - alpha = x(m+1:2*m); - - picture( param, FUNstr, BASISstr, BOUNDARYstr, ... - x, error_norm, t, alpha ); - - end; % if - - [param] = stopmenu(param, relative_error ); - - end; % while - - % extract various variables from x - [length_t, t, alpha, r, lower_bound, lambda] =... - def_var(param, x); - - disp( [ 'End COCA' ] ); -return; //GO.SYSIN DD coca.m echo compnorm.m 1>&2 sed >compnorm.m <<'//GO.SYSIN DD compnorm.m' 's/^-//' - function [error_norm, relative_error, t_max, alpha_max ] = ... - compnorm( param, FUNstr, BASISstr, BOUNDARYstr, x ); -%function [error_norm, relative_error, t_max, alpha_max ] = ... -%compnorm( param, FUNstr, BASISstr, BOUNDARYstr, x ); -% -% COMPNORM.M 02.07.92 -% -% Finds the maximum of the error function, starting from -% a discrete grid (defined by param(5)). -% -% Single_exchange = param(4) = 1 -% Returns the extremal point (t_max, alpha_max) that maximizes -% the error function (see: ERRFUN.M), the maximum of the -% error function, and the relative distance between this maximum -% and the lower bound. -% -% multiple_exchange = param(4) ~= 1 -% Returns the points (t_max, alpha_max) which have function values -% sufficiently close to the maximum of the error function -% (see: ERRFUN.M), the maximum of the error function, and -% the relative distance between this maximum and the lower bound. -% -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - C_dim = param(1); - R_dim = param(2); - single_exchange = param(4); - step = 1/param(5); - output = param(17); - -% extract various variables from x - [length_t, t, alpha, r, lower_bound, lambda] = ... - def_var( param, x ); - -% set grid defined by param(5) - dt = (0:step:1)'; - -% add old points - for j=1:length_t - I = find(t(j) == dt ); - if length(I) == 0 - dt = [dt;t(j)]; - end; % if - end; % for - -% add critical-points - if length( param ) > 20 - critical_points = param(21:length(param)); - for j=1:length(critical_points) - I=find(critical_points(j)==dt); - if length(I) == 0 - dt=[dt;critical_points(j)]; - end; % if - end; % for - end; % if - -% add end points - I = find(0 == dt); - if length(I) == 0 dt = [dt;0]; end; - I = find(1 == dt); - if length(I) == 0 dt = [dt;1]; end; - -% evaluate error function - dt = sort( dt ); - [error abs_error] = errfun( FUNstr, BASISstr, BOUNDARYstr, ... - C_dim, R_dim, lambda, dt ); - -% find (continuous) maximum of abs(error) - - % locate local (discrete) maximum of abs(error) - signum = sign_spe( abs_error(2:length(dt)) ... - -abs_error(1:length(dt)-1) ); - signum = signum(1:length(signum)-1) ... - + 2 * signum(2:length(signum)); - I = find( signum == -1 ); - - % define data for parabolic interpolation - abs_minus = abs_error( I ); - abs_middle = abs_error( I + 1 ); - abs_plus = abs_error( I + 2 ); - t_minus = dt( I ); - t_middle = dt( I + 1 ); - t_plus = dt( I + 2 ); - - % interpolate - d1 = ( abs_middle - abs_minus )./( t_middle - t_minus ); - d2 = ( abs_middle - abs_plus )./( t_middle - t_plus ); - d3 = ( d1 - d2 )./( t_minus - t_plus ); - % compute maximum of the parabola - t_int = ( t_minus + t_middle - d1./d3 ) / 2; - - % add maxima of the interpolation process - if length( t_int ) > 0 - [int_error, abs_int_error] =... - errfun( FUNstr, BASISstr, BOUNDARYstr, ... - C_dim, R_dim, lambda, t_int ); - - dt = [dt; t_int]; - error = [error; int_error]; - abs_error = [abs_error; abs_int_error]; - end; % if - - [error_norm,j] = max( abs_error ); - relative_error = ( error_norm - lower_bound ) / abs ( lower_bound ); - - if single_exchange == 1 - - error_max = error_norm; - t_max = dt(j); - alpha_max = angle( error(j) ); - - else % multiple exchange - - gap = error_norm - lower_bound; - I = find( abs_error >= (error_norm - gap/3 )); - t_max = dt( I ); - error = error( I ); - error_max = abs_error( I ); - [error_max,I] = sort(error_max); - t_max = t_max(I); - error = error( I ); - alpha_max = angle( error ); - - end % if - - - if output > 2 & single_exchange == 1 - disp( [ 't_max = ', num2str(t_max) ] ); - disp( [ 'alpha_max = ', num2str(alpha_max) ] ); - end; - - if output > 0 - disp( [ 'error_norm = ', num2str(error_norm) ] ); - disp( [ 'low. bound = ', num2str(lower_bound) ] ); - disp( [ 'rel. error = ', num2str(relative_error) ] ); - end - - if output > 3 - pause - end -return; //GO.SYSIN DD compnorm.m echo convstr.m 1>&2 sed >convstr.m <<'//GO.SYSIN DD convstr.m' 's/^-//' - function [ FUNstr, BASISstr, BOUNDARYstr ]=convstr( FUN, FUNarg,... - BASISfun, BASISarg, BOUNDARY, BOUNDARYarg, out ); -%function [ FUNstr, BASISstr, BOUNDARYstr ]=convstr( FUN, FUNarg,... -% BASISfun, BASISarg, BOUNDARY, BOUNDARYarg, out ); -% -% CONVSTR.M 16.07.92 -% -% Appends to the strings FUNstr, BASISstr, and BOUNDARYstr -% their default arguments (see below) and possible additional -% arguments defined in FUNarg, BASISarg, and BOUNDARYarg, respectively. -% -% default input and output arguments: -% [gamma, dgamma] = BOUNDARY(t) -% [f, df] = FUN(gamma) -% [phi, dphi] = BASISfun(gamma, C_dim) -% -% Note: in this version complex input parameters have to be -% splitted into their real and imaginary part. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% FUN - FUNstr = [FUN, '(gamma']; - for i=1:length(FUNarg) - eval(['P=FUNarg(',num2str(i),');']); - str=sprintf('%.16g',P); - FUNstr = [FUNstr,',',str]; - end; % for - FUNstr = [FUNstr, ')']; - -% BASIS - if length(BASISarg) == 0 - BASISstr = [BASISfun,'(gamma,C_dim)']; - else - BASISstr = [BASISfun,'(gamma,C_dim,[']; - for i=1:length(BASISarg) - eval(['P=BASISarg(',num2str(i),');']); - str=sprintf('%.16g',P); - if i==1 - BASISstr = [BASISstr,str]; - elseif i > 1 - BASISstr = [BASISstr,',',str]; - end % if - end % for - BASISstr = [BASISstr, '])']; - end % if - -% BOUNDARY - BOUNDARYstr = [BOUNDARY, '(t']; - for i=1:length(BOUNDARYarg) - eval(['P=BOUNDARYarg(',num2str(i),');']); - str=sprintf('%.16g',P); - BOUNDARYstr = [BOUNDARYstr,',',str]; - end % for - BOUNDARYstr = [BOUNDARYstr, ')']; - - if out > 1 - FUNstr - BASISstr - BOUNDARYstr - if out > 3 - pause - end; % if - end; % if -return; //GO.SYSIN DD convstr.m echo def_var.m 1>&2 sed >def_var.m <<'//GO.SYSIN DD def_var.m' 's/^-//' - function [length_t, t, alpha, r, lower_bound, lambda] =... - def_var( param, x ); -%function [length_t, t, alpha, r, lower_bound, lambda] =... -%def_var( param, x ); -% -% DEF_VAR.M 16.07.92 -% -% Extracts variables from x for usage in various functions. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - C_dim = param(1); - R_dim = param(2); - l = length(x); - length_t = (l - (R_dim + 1)) / 3; - t = x( 1: length_t); - alpha = x( length_t+1:2*length_t); - r = x(2*length_t+1:3*length_t); - lower_bound = x(3*length_t+1); - lambda = x(3*length_t+2:l); -return; //GO.SYSIN DD def_var.m echo ellips_c.m 1>&2 sed >ellips_c.m <<'//GO.SYSIN DD ellips_c.m' 's/^-//' - function [gamma, dgamma] = ellips_c( t, R, Rec, Imc ); -%function [gamma, dgamma] = ellips_c( t, R, Rec, Imc ); -% -% ELLIPS_C.M 16.07.92 -% -% Chebyshev ellipse with centre (Rec + i*Imc). -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - t = 2*pi*t; - gamma = ( (R+1/R) * cos(t) + i*(R-1/R) * sin(t) ) / 2 ... - + Rec + i * Imc; - dgamma = (-(R+1/R) * sin(t) + i*(R-1/R) * cos(t) ) * pi; -return; //GO.SYSIN DD ellips_c.m echo ellipse.m 1>&2 sed >ellipse.m <<'//GO.SYSIN DD ellipse.m' 's/^-//' - function [gamma, dgamma] = ellipse( t, R ); -%function [gamma, dgamma] = ellipse( t, R ); -% -% ELLIPSE.M 16.07.92 -% -% Chebyshev ellipse. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - t = 2*pi*t; - gamma = ( (R+1/R) * cos(t) + i*(R-1/R) * sin(t) ) / 2; - dgamma = (-(R+1/R) * sin(t) + i*(R-1/R) * cos(t) ) * pi; -return; //GO.SYSIN DD ellipse.m echo errfun.m 1>&2 sed >errfun.m <<'//GO.SYSIN DD errfun.m' 's/^-//' - function [error, abs_error] = ... - errfun( FUNstr, BASISstr, BOUNDARYstr, C_dim, R_dim, lambda, t ) -%function [error, abs_error] = ... -%errfun( FUNstr, BASISstr, BOUNDARYstr, C_dim, R_dim, lambda, t ) -% -% ERRFUN.M 16.07.92 -% -% Computes the error-function: -% error = f - phi.' * lambda, -% abs_error = abs( error ). -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. - - [gamma,dgamma] = eval( BOUNDARYstr ); - [phi,dphi] = eval( BASISstr ); - -% note, for complex coefficients -% phi(C_dim + j) = i * phi(j) - - if R_dim ~= C_dim - phi = [ phi; i * phi]; - dphi = [dphi; i * dphi]; - end - - [f,df] = eval( FUNstr ); - error = f - phi.' * lambda; - abs_error = abs( error ); - -return; - - //GO.SYSIN DD errfun.m echo exam_11.m 1>&2 sed >exam_11.m <<'//GO.SYSIN DD exam_11.m' 's/^-//' -% EXAM_11.M 17.07.92 -% -% Approximate f(z) = cos(z), -% by phi(j,z) = z^j, j=0,2,4,6, -% on the unit circle. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - clear - format short - -% name of approximant - FUNstr = 'fcos'; - -% this function has no additional parameter - FUNarg = []; - -% name of basis function - BASISstr = 'monom'; - -% this is the special choice - BASISarg = [0 2 4 6]; - -% name of boundary - BOUNDARYstr = 'circle'; - -% this boundary needs no parameter - BOUNDARYarg = []; - -% the function MONOM.M has no additional parameter, so the -% complex dimension can be defined by the length of the arguments - C_dim = length(BASISarg); - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter -% for further information on parameters see SET_PARA.M - param(17) = 4; % for very extended output and pauses - param(18) = 2; % plots are shown at each iteration step - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_11.m echo exam_12.m 1>&2 sed >exam_12.m <<'//GO.SYSIN DD exam_12.m' 's/^-//' -% EXAM_12.M 17.07.92 -% -% Approximate f(z) = sin(z), -% by phi(j,z) = z^j, j=1,2,3, -% on the unit circle. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - -% name of approximant - FUNstr = 'fsin'; - -% this function has no additional parameter - FUNarg = []; - -% name of basis function - BASISstr = 'monom'; - -% this is the special choice - BASISarg = 1:3; - -% name of boundary - BOUNDARYstr = 'circle'; - -% this boundary needs no parameter - BOUNDARYarg = []; - -% the function MONOM.M has no additional parameter, so the -% complex dimension can be defined by the length of the arguments - C_dim = length(BASISarg); - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are complex - real_coef = 0; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter -% for further information on parameters see SET_PARA.M - param(17) = 4; % for very extended output and pauses - param(18) = 2; % plots are shown at each iteration step - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_12.m echo exam_13.m 1>&2 sed >exam_13.m <<'//GO.SYSIN DD exam_13.m' 's/^-//' -% EXAM_13.M 17.07.92 -% -% Approximate f(z) = exp(z), -% by phi(j,z) = z^j, j=0,1,2,3, -% on the unit quadrat. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - -% name of approximant - FUNstr = 'fexp'; - -% this function has no additional parameter - FUNarg = []; - -% name of basis function - BASISstr = 'monom'; - -% this is the special choice - BASISarg = [0 1 2 3]; - -% name of boundary - BOUNDARYstr = 'recangle'; - -% this boundary needs the parameters (see: RECANGLE.M) - right = 1; % right coordinate of quadrat - upper = 1; % upper coordinate of quadrat - left = -1; % left coordinate of quadrat - lower = -1; % lower coordinate of quadrat - BOUNDARYarg = [right upper left lower]; - -% the function MONOM.M has no additional parameter, so the -% complex dimension can be defined by the length of the arguments - C_dim = length(BASISarg); - -% here we define critical points, i.e. corners -% (the partial derivative does not exists) - critical_points = [0, 0.25, 0.5, 0.75, 1]'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter -% for further information on parameters see SET_PARA.M - param(17) = 4; % for very extended output and pauses - param(18) = 2; % plots are shown at each iteration step - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_13.m echo exam_14.m 1>&2 sed >exam_14.m <<'//GO.SYSIN DD exam_14.m' 's/^-//' -% EXAM_14.M 17.07.92 -% -% Approximate f(z) = z^5, -% by phi(j,z) = z^j, j=0:4, -% on the real unit interval. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - clear - format short - -% name of approximant - FUNstr = 'fmonom'; - -% this function has no additional parameter - FUNarg = [5]; - -% name of basis function - BASISstr = 'monom'; - -% this is the special choice - BASISarg = [0:4]; - -% name of boundary - BOUNDARYstr = 'interval'; - -% this boundary needs no parameter - BOUNDARYarg = []; - -% the function MONOM.M has no additional parameter, so the -% complex dimension can be defined by the length of the arguments - C_dim = length(BASISarg); - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter -% for further information on parameters see SET_PARA.M - param(4) = 0; % multiple exchange in IMPROVE.M - param(17) = 4; % for very extended output and pauses - param(18) = 2; % plots are shown at each iteration step - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_14.m echo exam_21.m 1>&2 sed >exam_21.m <<'//GO.SYSIN DD exam_21.m' 's/^-//' -% EXAM_21.M 17.07.92 -% -% Approximate f(z) = 1, -% by phi(j,z) = (z-c)^j, j=1,2,3,4,5, -% where c is defined as below, -% -% on a chebyshev ellipse with foci 1, -1, -% semi axis (R + 1/R)/2, (R - 1/R)/2; -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - - R = 3; - Rr = 3.001; - Rec = (Rr + 1/Rr)/2; - Imc = 0; - % shift parameter for - % basis function: c = Rec + i*Imc; - -% name of approximant - FUNstr = 'fone'; - -% this function has no additional parameter - FUNarg = []; - -% name of basis function - BASISstr = 'monom_c'; - -% this is the special choice - BASISarg = [1:5, Rec, Imc]; - -% name of boundary - BOUNDARYstr = 'ellipse'; - -% the boundary parameter - BOUNDARYarg = [R]; - -% the function MONOM_C.M has two additional parameters, so the -% complex dimension must be defined directly - C_dim = 5; - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 3; % for very extended output - param(18) = 1; % plots are shown after remez iteration - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_21.m echo exam_22.m 1>&2 sed >exam_22.m <<'//GO.SYSIN DD exam_22.m' 's/^-//' -% EXAM_22.M 17.07.92 -% -% Approximate f(z) = 1, -% by phi(j,z) = z^j, j=1,2,3,4,5, -% -% on a chebyshev ellipse with foci 1, -1, -% semi axis (R + 1/R)/2, (R - 1/R)/2, and centre rec + i*Imc; -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - - R = 3; - Rec = 1.67; - Imc = 0; - -% name of approximant - FUNstr = 'fone'; - -% this function has no additional parameter - FUNarg = []; - -% name of basis function - BASISstr = 'monom'; - -% this is the special choice - BASISarg = [1:5]; - -% name of boundary - BOUNDARYstr = 'ellips_c'; - -% the boundary parameter - BOUNDARYarg = [R, Rec, Imc]; - -% the function MONOM.M has no additional parameter, so the -% complex dimension can be defined by the length of the arguments - C_dim = length(BASISarg); - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 3; % for very extended output - param(18) = 1; % plots are shown after remez iteration - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_22.m echo exam_23.m 1>&2 sed >exam_23.m <<'//GO.SYSIN DD exam_23.m' 's/^-//' -% EXAM_22.M 17.07.92 -% -% Approximate f(z) = z^4, -% by phi(j,z) = z^j, j=0,1,2,3, -% on the unit disk. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - -% name of approximant - FUNstr = 'fmonom'; - -% this function has one additional parameter - FUNarg = [4]; - -% name of basis function - BASISstr = 'monom'; - -% this is the special choice - BASISarg = [0:3]; - -% name of boundary - BOUNDARYstr = 'circle'; - -% the boundary parameter - BOUNDARYarg = []; - -% the function MONOM.M has no additional parameter, so the -% complex dimension can be defined by the length of the arguments - C_dim = length(BASISarg); - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 3; % for very extended output - param(18) = 1; % plots are shown after remez iteration - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_23.m echo exam_31.m 1>&2 sed >exam_31.m <<'//GO.SYSIN DD exam_31.m' 's/^-//' -% EXAM_31.M 17.07.92 -% -% Approximate f(z) = 1/(z-c), -% where c is defined below, -% by phi(j,z) = z^j, j=0,1,2, -% on the unit circle. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - clear - format short - -% define parameter c = Rec + i*Imc; - Rec = 2; - Imc = 1; - -% name of approximant and additional parameter - FUNstr = 'fpole_c'; - FUNarg = [Rec, Imc]; - -% name of basis function and parameter - BASISstr = 'monom'; - BASISarg = [0, 1, 2]; - -% name of boundary and parameter - BOUNDARYstr = 'circle'; - BOUNDARYarg = []; - - C_dim = 3; - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are complex - real_coef = 0; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 2; % for middle output - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_31.m echo exam_32.m 1>&2 sed >exam_32.m <<'//GO.SYSIN DD exam_32.m' 's/^-//' -% EXAM_32.M 17.07.92 -% -% Approximate f(z) =z^5, -% by phi(j,z) = z^j, j=1,3, -% on a circular sector (angle 30 degree). -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - clear - format short - -% define parameter c = Rec + i*Imc; - Rec = 2; - Imc = 1; - -% name of approximant and additional parameter - FUNstr = 'fmonom'; - FUNarg = [5]; - -% name of basis function and parameter - BASISstr = 'monom'; - BASISarg = [1,3]; - -% name of boundary and parameter - BOUNDARYstr = 'sector'; - BOUNDARYarg = [30]; - - C_dim = length( BASISarg ); - -% here we define critical points, i.e. corners -% (the partial derivative does not exists) - critical_points = [0, 0.5, 1]'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 2; % for middle output - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_32.m echo exam_33.m 1>&2 sed >exam_33.m <<'//GO.SYSIN DD exam_33.m' 's/^-//' -% EXAM_33.M 17.07.92 -% -% Approximate f(z) =z^5, -% by phi(j,z) = z^j, j=1,3, -% on a rectangle with corners defined below. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - clear - format short - -% define parameter c = Rec + i*Imc; - Rec = 2; - Imc = 1; - -% name of approximant and additional parameter - FUNstr = 'fmonom'; - FUNarg = [5]; - -% name of basis function and parameter - BASISstr = 'monom'; - BASISarg = [1,3]; - -% name of boundary and parameter - BOUNDARYstr = 'recangle'; - right = 2; % right coordinate of quadrat - upper = 1; % upper coordinate of quadrat - left = -2; % left coordinate of quadrat - lower = -1; % lower coordinate of quadrat - BOUNDARYarg = [right upper left lower]; - - C_dim = length( BASISarg ); - -% here we define critical points, i.e. corners -% (the partial derivative does not exists) - critical_points = [0, 0.25, 0.5, 0.75, 1]'; - -% we will assume, that the coefficients are real - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 2; % for middle output - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_33.m echo exam_34.m 1>&2 sed >exam_34.m <<'//GO.SYSIN DD exam_34.m' 's/^-//' -% EXAM_34.M 17.07.92 -% -% Approximate f(z) = 1/(z-c), -% where c is defined below, -% by phi(j,z) = z^j, j=0,1,2, -% on a circle shifted by cc = Recc + i*Imcc with radius r = 0.5. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - clear - format short - -% define parameters c, cc, r - Rec = 2; - Imc = 1; - Recc = -1; - Imcc = -1; - r = 0.5; - -% name of approximant and additional parameter - FUNstr = 'fpole_c'; - FUNarg = [Rec, Imc]; - -% name of basis function and parameter - BASISstr = 'monom'; - BASISarg = [0, 1, 2]; - -% name of boundary and parameter - BOUNDARYstr = 'circle_c'; - BOUNDARYarg = [Recc, Imcc, r]; - - C_dim = 3; - -% there are no points, at which the partial derivative does not exists - critical_points = []'; - -% we will assume, that the coefficients are complex - real_coef = 0; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 2; % for middle output - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_34.m echo exam_41.m 1>&2 sed >exam_41.m <<'//GO.SYSIN DD exam_41.m' 's/^-//' -% EXAM_41.M 17.07.92 -% -% Approximate f(z) = 1, -% by phi(j,z) = z^j, j=1:10, -% -% on a chebyshev ellipse with foci 1, -1, -% semi axis (R + 1/R)/2, (R - 1/R)/2, and -% centre Rec + i*Imc, parameters see below. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - - clear - format short - - R = 3; - Rec = 1.67; - Imc = 0; - - FUNstr = 'fone'; - FUNarg = []; - BASISstr = 'monom'; - BASISarg = [1:10]; - BOUNDARYstr = 'ellips_c'; - BOUNDARYarg = [R, Rec, Imc]; - C_dim = length(BASISarg); - critical_points = []'; - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_41.m echo exam_42.m 1>&2 sed >exam_42.m <<'//GO.SYSIN DD exam_42.m' 's/^-//' -% EXAM_42.M 17.07.92 -% -% Approximate f(z) = 1/(z-c), -% where c is defined below, -% by phi(j,z) = z^j, j=0:5, -% on the unit circle. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - - clear - format short - - Rec = 2; - Imc = 1; - - FUNstr = 'fpole_c'; - FUNarg = [Rec Imc]; - BASISstr = 'monom'; - BASISarg = [0:5]; - BOUNDARYstr = 'circle'; - BOUNDARYarg = []; - C_dim = length(BASISarg); - critical_points = []'; - real_coef = 0; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_42.m echo exam_43.m 1>&2 sed >exam_43.m <<'//GO.SYSIN DD exam_43.m' 's/^-//' -% EXAM_43.M 17.07.92 -% -% Approximate f(z) = z, -% by phi(j,z) = z^j, j=2:10, -% on a L-shaped region, -% see Example 2 out of Glashoff/Roleff, Math. Comp. 36:233-239, 1981. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - - d = 2/sqrt(5); - FUNstr = 'fmonom'; - FUNarg = [1]; - BASISstr = 'monom'; - BASISarg = 2:9; - BOUNDARYstr = 'l_shape'; - BOUNDARYarg = [d]; - C_dim = length(BASISarg); - critical_points = [0, 0.25, 0.5, 1]'; - real_coef = 1; - -% define the default parameter - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameters - param(5) = 100; % number of steps in COMPNORM.M - param(6) = 1e-2; % bound for relative error -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_43.m echo exam_51.m 1>&2 sed >exam_51.m <<'//GO.SYSIN DD exam_51.m' 's/^-//' -% EXAM_51.M 17.07.92 -% -% Approximate f(z) = z, -% by phi(j,z) = z^j, j=2:10, -% -% on a chebyshev ellipse with foci 1, -1, -% semi axis (R + 1/R)/2, (R - 1/R)/2, and -% centre Rec + i*Imc, parameters see below. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - - clear - format short - - R = 3; - Rec = 1.67; - Imc = 0; - - FUNstr = 'fmonom'; - FUNarg = [1]; - BASISstr = 'monom'; - BASISarg = [2:10]; - BOUNDARYstr = 'ellips_c'; - BOUNDARYarg = [R, Rec, Imc]; - C_dim = length(BASISarg); - critical_points = []'; - real_coef = 1; - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(17) = 0; % no output - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_51.m echo exam_52.m 1>&2 sed >exam_52.m <<'//GO.SYSIN DD exam_52.m' 's/^-//' -% EXAM_52.M 17.07.92 -% -% Approximate f(z) = 1, -% by phi(j,z) = z^j, j=1:10, -% -% on two disjoint cirles, with centres Rem(i), Imm(i) and -% radius r(i), i=1,2. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - - Rem1 = -2; Imm1 = -1; r1 = 1;, - Rem2 = 1; Imm2 = 0.5; r2 = 0.5; - - FUNstr = 'fmonom'; - FUNarg = [1]; - BASISstr = 'monom'; - BASISarg = [2:10]; - BOUNDARYstr = 'two_circ'; - BOUNDARYarg = [Rem1, Imm1, r1, Rem2, Imm2, r2]; - C_dim = length(BASISarg); - critical_points = []'; - real_coef = 0; - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(5) = 100; % numer of steps in COMPNORM.M - param(6) = 1e-2; % relative error bound - param(17) = 3; % extended output - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_52.m echo exam_53.m 1>&2 sed >exam_53.m <<'//GO.SYSIN DD exam_53.m' 's/^-//' -% EXAM_53.M 17.07.92 -% -% Approximate f(z) = exp(z), -% by phi(j,z) = z^j, j=0:5, -% -% on the unit circle. -% -% This program is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - clear - format short - - FUNstr = 'fexp'; - FUNarg = []; - BASISstr = 'monom'; - BASISarg = [0:5]; - BOUNDARYstr = 'circle'; - BOUNDARYarg = []; - C_dim = length(BASISarg); - critical_points = []'; - real_coef = 1; - [param] = set_para( critical_points, C_dim, real_coef ); - -% define special parameter - param(4) = 0; % multiple exchange in IMPROVE.M - param(6) = 1e-2; % relative error bound - param(17) = 3; % extended output - param(18) = 2; % plots after each iteration step - -% call the COCA package - [param, t, alpha, r, lower_bound, lambda, error_norm ] =... - coca(param, FUNstr, FUNarg, BASISstr, BASISarg,... - BOUNDARYstr, BOUNDARYarg); - -% some output - disp( [ ' t', ' alpha',' r'] ) - disp([t, alpha, r]); - disp(['lambda']); - disp([lambda]); //GO.SYSIN DD exam_53.m echo fcos.m 1>&2 sed >fcos.m <<'//GO.SYSIN DD fcos.m' 's/^-//' - function [f, df] = fcos( gamma ); -%function [f, df] = fcos( gamma ); -% -% FCOS.M 16.07.92 -% -% f(z) = cos(z). -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - f = cos(gamma); - df = -sin(gamma); -return; //GO.SYSIN DD fcos.m echo fexp.m 1>&2 sed >fexp.m <<'//GO.SYSIN DD fexp.m' 's/^-//' - function [f, df] = fexp( gamma ); -%function [f, df] = fexp( gamma ); -% -% FEXP.M 16.07.92 -% -% f(z) = exp(z). -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - f = exp(gamma); - df = -exp(gamma); -return; //GO.SYSIN DD fexp.m echo fmonom.m 1>&2 sed >fmonom.m <<'//GO.SYSIN DD fmonom.m' 's/^-//' - function [f, df] = fmonom( gamma , N ); -%function [f, df] = fmonom( gamma , N ); -% -% FMONOM.M 15.07.92 -% -% f(z) = z^N. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - f = gamma.^N; - df = N*gamma.^(N-1); -return; - - - - //GO.SYSIN DD fmonom.m echo fone.m 1>&2 sed >fone.m <<'//GO.SYSIN DD fone.m' 's/^-//' - function [f, df] = fone( gamma ); -%function [f, df] = fone( gamma ); -% -% FONE.M 16.07.92 -% -% f(z) = 1. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - f = ones(gamma); - df = zeros(gamma); -return; //GO.SYSIN DD fone.m echo fpole_c.m 1>&2 sed >fpole_c.m <<'//GO.SYSIN DD fpole_c.m' 's/^-//' - function [f, df] = fpole_c( z, Rec, Imc ); -%function [f, df] = fpole_c( z, Rec, Imc ); -% -% FPOLE_C.M 16.07.92 -% -% 1 -% f(z) = -----, c = Rec + i * Imc. -% z - c -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - c = Rec + sqrt(-1) * Imc; - f = ones(length(z), 1)./(z - c); - df = -ones(length(z), 1)./(z - c).^2; -return; //GO.SYSIN DD fpole_c.m echo fsin.m 1>&2 sed >fsin.m <<'//GO.SYSIN DD fsin.m' 's/^-//' - function [f, df] = fsin( gamma ); -%function [f, df] = fsin( gamma ); -% -% FSIN.M 16.07.92 -% -% f(z) = sin(z). -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - f = sin(gamma); - df = cos(gamma); -return; //GO.SYSIN DD fsin.m echo improve.m 1>&2 sed >improve.m <<'//GO.SYSIN DD improve.m' 's/^-//' - function [x, A, c] = ... - improve(param, FUNstr, BASISstr, BOUNDARYstr, ... - x, t_max, alpha_max, A, c); -%function [x, A, c] = ... -%improve(param, FUNstr, BASISstr, BOUNDARYstr, ... -% x, t_max, alpha_max, A, c); -% -% IMPROVE.M 16.07.92 -% -% Improves the lower bound for the minimal deviation, by -% solving the dual formulation of the (discrete) problem with -% the simplex method. -% -% Single_exchange = param(4) = 1 -% One simplex step is performed directly. -% -% multiple_exchange = param(4) ~= 1 -% Call of SIMPLEX.M. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - C_dim = param(1); - R_dim = param(2); - single_exchange = param(4); - output = param(17); - critical_points = param(20:length(param)); - -% extract various variables from x - [m, t, alpha, r, lower_bound, lambda] =... - def_var( param, x ); - - y = [lower_bound; lambda]; - - if output > 1 - disp([' ']) - disp(['IMPROVE']) - end - - if single_exchange == 1 - - % one simplex step 'by hand' - - % compute new column with respect to (t_max,alpha_max) - [new_a, new_c] = ... - set_colm( param, FUNstr, BASISstr, BOUNDARYstr, ... - t_max, alpha_max ); - - % find pivot element p - % zeta = r(p) / d(p) = min{ r(i) / d(i), d(i) > 0 } - d = A\new_a; - I = find( d > 0 ); - [dummy,p] = min( r(I)./d(I) ); - p = I(p); - - % exchange - t(p) = t_max; - alpha(p) = alpha_max; - - % new matrix A and right hand side c - [A,c] = ... - set_colm( param, FUNstr, BASISstr, BOUNDARYstr, t, alpha ); - - % new weights - b = eye( m, 1 ); - r = A\b; - - else % multiple_exchange - - t = [t; t_max ]; - alpha = [alpha; alpha_max ]; - r = [r; zeros(t_max)]; - - % define matrix for SIMPLEX - [A,c] = set_colm( param, FUNstr, BASISstr, BOUNDARYstr,... - t, alpha ); - b = eye( m, 1 ); - - % call SIMPLEX - [r,I] = simplex( -c, A, b, 1:m ); - - % extract active values - t = t(I); - alpha = alpha(I); - A = A(:,I); - r = r(I); - c = c(I); - - end; %if - - if output > 2 - disp( [ ' t', ' alpha', ... - ' r', ' lambda' ] ) - disp([t, alpha, r, y ]); - if output > 3 - pause; - end; %if - end; %if - -% collect x - x = [t; alpha; r; y]; -return; //GO.SYSIN DD improve.m echo initial.m 1>&2 sed >initial.m <<'//GO.SYSIN DD initial.m' 's/^-//' - function [ x, A, c ] = initial( param, FUNstr, BASISstr, BOUNDARYstr, x ); -%function [ x, A, c ] = initial( param, FUNstr, BASISstr, BOUNDARYstr, x ); -% -% INITIAL.M 16.07.92 -% -% Computes the initial extremal points (see: SET_EXTR.M), -% matrix A, right hand side c, and weights r. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - C_dim = param(1); - R_dim = param(2); - initial_extremal_points = param(3); - output = param(17); - m = R_dim + 1; - t = x( 1: m); - alpha = x(m+1:2*m); - - if output > 1 - disp([' ']) - disp(['INITIALIZATION']) - end - - if output > 2 - disp(['Complex-Dimension = ', int2str( C_dim )]); - disp(['Real -Dimension = ', int2str( R_dim )]); - end; - -% set initial extremal points - if initial_extremal_points ~= 0 - [t, alpha] = set_extr( param, x ); - end; - -% set system matrix A and right hand side c - [A,c] = set_colm( param, FUNstr, BASISstr, BOUNDARYstr, t, alpha ); - -% A numerically singular ? - if min(svd(A)) < 2*eps - disp(['A is numerically singular']); - error('try different starting values (good luck)') - end; - -% compute weights - b = eye( m, 1 ); - r = A\b; - -% negative initial weights ? -% if weigth r(i) negative, choose angle alpha(i) = alpha(i) + pi - I = find( r < 0 ); - if length(I) > 0 - - alpha(I) = alpha(I) + pi; - I = find( alpha >= pi ); - if length(I) > 0 - alpha(I) = alpha(I) - 2 * pi; - end; % if - - % note: the following redefinition of A - % and c is not necessary: - % only columns related to modified alpha(i)'s - % have to be multiplicated by -1 - % (except for first row which still remains 1). - - [A,c] = set_colm( param, FUNstr, BASISstr, BOUNDARYstr, t, alpha ); - - % An update technique can be used instead of - % solving the linear system A * r = b - - r = A\b; - - end % if - - if output > 2 - disp([' t',' alpha',' c',' r']) - disp([t alpha c r]) - if output > 3 - pause - end; - end; %if - -% collect x - x( 1: m) = t; - x( m+1:2*m) = alpha; - x(2*m+1:3*m) = r; -return; //GO.SYSIN DD initial.m echo interval.m 1>&2 sed >interval.m <<'//GO.SYSIN DD interval.m' 's/^-//' - function [gamma, dgamma] = interval( t ); -%function [gamma, dgamma] = interval( t ); -% -% INTERVAL.M 16.07.92 -% -% Unit-interval. -% t in [0,1] -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - t = 2*t-1; - gamma = t; - dgamma = 2*ones(t); -return; //GO.SYSIN DD interval.m echo l_shape.m 1>&2 sed >l_shape.m <<'//GO.SYSIN DD l_shape.m' 's/^-//' - function [gamma, dgamma] = l_shape( t, d ); -%function [gamma, dgamma] = l_shape( t, d ); -% -% L_SHAPE.M 16.07.92 -% -% Example 2 out of Glashoff/Roleff, Math. Comp. 36:233-239, 1981. -% Note: t in [0,1] in our example. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - [dummy, col] = size(t); - if col ~= 1 - t = t'; - end; - - t = 2 * t; - gamma = zeros(length(t),1); - dgamma = ones(length(t),1); - i = sqrt(-1); - - I1 = find( t <= 0.5 ); - if length(I1) > 0 - gamma(I1) = 1/2 + t(I1) + i*t(I1); - dgamma(I1) = 2*ones(length(I1),1) + i*2*ones(length(I1),1); - end; % if - - I2 = find( t > 0.5 & t <= 1 ); - if length(I2) > 0 - gamma(I2) = 3/2 - t(I2) + i * t(I2); - dgamma(I2) = -2*ones(length(I2),1) + i*2*ones(length(I2),1); - end; - - I3 = find( t > 1 & t <= 2 ); - if length(I3) > 0 - gamma(I3) = 3/2 - t(I3) + i*(2-t(I3)); - dgamma(I3) = -2*ones(length(I3),1) - i*2*ones(length(I3),1); - end; - - gamma = d * gamma; - dgamma = d * dgamma; - -return; //GO.SYSIN DD l_shape.m echo menustr.m 1>&2 sed >menustr.m <<'//GO.SYSIN DD menustr.m' 's/^-//' -function kstr = ... -menu( s0, s1, s2, s3, s4, ... - s5, s6, s7, s8, s9, ... - s10, s11, s12, s13, s14 ); -% -% MENUSTR.M 05.07.92 -% -% Note: this routine is basically the program MENU.M out -% of the Matlab toolbox, but it returns a string instead -% of a number. -% -% This function is part of the COCA - package. -% - disp(' ') - disp(['----- ',s0,' -----']) - disp(' ') - for i=1:(nargin-1) - disp([' ',int2str(i),') ',eval(['s',int2str(i)])]) - end - disp(' ') - k = input( 'Select a menu number: ' ); - kstr = eval( ['s',int2str(k) ] ); -return; //GO.SYSIN DD menustr.m echo monom.m 1>&2 sed >monom.m <<'//GO.SYSIN DD monom.m' 's/^-//' - function [phi,dphi] = monom( z, C_dim, choice ); -%function [phi,dphi] = monom( z, C_dim, choice ); -% -% MONOM.M 15.07.92 -% -% Computes the (R_dim x length(z)) matrices phi, dphi, -% with i-th column -% -% phi(:,j) = z^ choice(j) -% dphi(:,j) = choice(j) * z^( choice(j) - 1 ) -% -% for complex or real vectors z=(z_1,z_2,...); -% choice is an arbitrary integer vector, e.g. [0 2 4]. -% -% For an additional shift parameter see MONOM_C.M -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - [row,dummy] = size(z); - if row > 1 - z = z.'; - end; - - m = length(z); - phi = zeros( C_dim, m); - dphi = zeros( C_dim, m); - h = ones( 1, m ); - dh = zeros( 1, m ); - j = 0; - k = 1; - - while k <= C_dim - if choice(k) == j - phi(k,:) = h; - dphi(k,:) = choice(k) * dh; - k = k + 1; - end; % if - dh = h; - h = z.*h; - j = j + 1; - - end; % while -return //GO.SYSIN DD monom.m echo monom_c.m 1>&2 sed >monom_c.m <<'//GO.SYSIN DD monom_c.m' 's/^-//' - function [phi,dphi] = monom_c( z, C_dim, choice ); -%function [phi,dphi] = monom_c( z, C_dim, choice ); -% -% MONOM_C.M 15.07.92 -% -% Computes the (R_dim x length(z)) matrices phi, dphi, -% with i-th column -% -% phi(:,j) = (z-c)^ choice(j) -% dphi(:,j) = choice(j) * (z-c)^( choice(j) - 1 ) -% -% for complex or real vectors z=(z_1,z_2,...); -% choice is a vector, defining the desired powers of (z-c) -% and the shift parameter c = Rec + i * Imc, e.g. [0 2 4 Rec Imc]. -% -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - [row,dummy] = size(z); - if row > 1 - z = z.'; - end; - - c = choice(length(choice) - 1) + i*choice(length(choice)); - m = length(z); - phi = zeros( C_dim, m); - dphi = zeros( C_dim, m); - h = ones( 1, m ); - dh = zeros( 1, m ); - j = 0; - k = 1; - while k <= C_dim - if choice(k) == j - phi(k,:) = h; - dphi(k,:) = choice(k) * dh; - k = k + 1; - end; % if - dh = h; - h = (z - c).*h; - j = j + 1; - - end; % while -return //GO.SYSIN DD monom_c.m echo newton_f.m 1>&2 sed >newton_f.m <<'//GO.SYSIN DD newton_f.m' 's/^-//' - function q = newton_f( param, FUNstr, BASISstr, BOUNDARYstr, x, length_t ); -%function q = newton_f( param, FUNstr, BASISstr, BOUNDARYstr, x, length_t ); -% -% NEWTON_F.M 16.07.92 -% -% This function defines the nonlinear function q for usage in -% Q_NEWTON.M. q consists of four parts: -% -% q1 = Real( error function * exp( -i*alpha ) ) - lower_bound -% q2 = Imag( error function * exp( -i*alpha ) ) -% q3 = Real( d/dt error function * exp( -i*alpha ) ) -% q4 = A( t(I), alpha(I) ) * r - b, I = 1:R_dim+1; -% -% -% If the partial derivative does not exist (q3) for the parameter value -% t(cp), i.e. dist( t(cp) - critical_points(j) ) is less then param(12), -% the equation -% Real( d/dt error function (t(cp)) * exp( -i*alpha(t(cp)) ) ) -% is substituted by -% t(cp) - critical_points(j) = 0 -% -% Note: this approach assumes, that all points with non existing -% derivative are collected in critical_points. -% -% For further information on the input parameter see REMEZ.M, SET_COLM.M, -% and SET_PARA.M. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. - - C_dim = param(1); - R_dim = param(2); - q = zeros(x,1); - - % extract various variables from x - [length_t, t, alpha, r, lower_bound, lambda] =... - def_var( param, x ); - - [gamma, dgamma] = eval( BOUNDARYstr ); - [phi,dphi] = eval( BASISstr ); - [f,df] = eval( FUNstr ); - - if R_dim ~= C_dim - phi = [phi; i * phi]; - dphi = [dphi; i * dphi]; - end - - error = f - phi.' * lambda; - derror = dgamma.*(df - dphi.' * lambda); - - - q( 1: length_t ) = ... - real( error.*exp(-i*alpha) ) - lower_bound; - q( length_t+1:2*length_t ) = imag( error.*exp(-i*alpha) ); - q( 2*length_t+1:3*length_t ) = real( derror.*exp(-i*alpha) ); - q( 3*length_t+1:length(x)-1) = real( phi*diag(exp(-i*alpha)) ) * r; - q(length(x)) = sum(r) - 1; - - [dummy,col] = size(q); - if col ~= 1 - q = q'; - end; - -% add critical-points - if length( param ) > 20 - - critical_points = param(21:length(param)); - - for j = 1:length_t - - k = find( abs( t(j) - critical_points ) < param(12) ); - if length(k) > 0 - if length(k) > 1 - error('param(12) too large !!'); - end % if - q(2*length_t+j) = t(j) - critical_points(k); - end % if - - end %for - - end % if -return; //GO.SYSIN DD newton_f.m echo picture.m 1>&2 sed >picture.m <<'//GO.SYSIN DD picture.m' 's/^-//' - function picture( param, FUNstr, BASISstr, BOUNDARYstr, x, ... - error_norm, t_max, alpha_max ); -%function picture( param, FUNstr, BASISstr, BOUNDARYstr, x, ... -% error_norm, t_max, alpha_max ); -% -% PICTURE.M 17.07.92 -% -% subplot(221): boundary, extremal points -% subplot(222): error function, extremal points, error norm circle -% subplot(223): abs(error function), extremal points, -% error norm, lower bound -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - C_dim = param(1); - R_dim = param(2); - critical_points = param(20:length(param)); - -% extract various variables from x - [length_t, t, alpha, r, lower_bound, lambda] =... - def_var( param, x ); - -% plot discretization - d = linspace( 0, 1, 500 )'; - -% corners should be ploted in any case - for j=1:length(critical_points) - I = find( critical_points(j) == d ); - if length(I) == 0, d = [d; critical_points(j)]; end - end %for - - d = sort(d); - [ error abs_error ] = ... - errfun( FUNstr, BASISstr, BOUNDARYstr, C_dim, R_dim, lambda, d ); - -% plot boundary - clg; - - dummy = t; - t = d; - gamma = eval( BOUNDARYstr ); - -% compute visually pleasing axis coordinates - maRe = max( real( gamma ) ); - miRe = min( real( gamma ) ); - maIm = max( imag( gamma ) ); - miIm = min( imag( gamma ) ); - - lRe = ( maRe - miRe ); - mRe = ( maRe + miRe )/2; - lIm = ( maIm - miIm ); - mIm = ( maIm + miIm )/2; - - if lRe > lIm - sg = 0.55 * lRe; - else - sg = 0.55 * lIm; - end %if - - -% plot boundary, extremal points - subplot(221); - axis( [mRe-sg, mRe+sg, mIm-sg, mIm+sg] ); - axis( 'square' ); - plot(real(gamma),imag(gamma),'w'); - hold; - t = dummy; - gamma = eval( BOUNDARYstr ); - plot( real(gamma),imag(gamma), 'ow' ); - plot([-100 100],[0 0],'w'); - plot([0 0],[-100 100],'w'); - hold - - -% plot error function, extremal points, error_norm circle - subplot(222) - axis('square'); - plot( error_norm * exp(i*d*2*pi),'w'), - hold - plot([-1000 1000],[0 0],'w'); - plot([0 0],[-1000 1000],'w'); - plot( real(error),imag(error), 'w' ); - plot( real( error_norm * exp(i*alpha_max)), ... - imag( error_norm * exp(i*alpha_max)) ,'ow') - hold - axis('normal'); - - -% plot abs(error function), extremal points, -% error norm, lower bound - Maxv = max( abs_error ); - minv = min( abs_error ); - MV = ( Maxv - minv )/20; - subplot(223) - axis([0 1 minv-MV Maxv+MV]); - plot( d, abs_error, 'w' ) - hold - plot( t_max, error_norm*ones(length(t_max)), 'ow' ) - plot([0 1],[lower_bound lower_bound],'w') - plot([0 1],[error_norm, error_norm],'w'); - for j=1:length(t) - plot([t(j) t(j)],[-1000 1000],'--w'); - end - xlabel(['lower b. = ',num2str(lower_bound),' ',... - 'norm = ',num2str(error_norm)]) ; - hold - - axis( 'normal' ); - pause; -return; //GO.SYSIN DD picture.m echo plot_ext.m 1>&2 sed >plot_ext.m <<'//GO.SYSIN DD plot_ext.m' 's/^-//' - function plot_ext( param, BOUNDARYstr, x, number ) -%function plot_ext( param, BOUNDARYstr, x, number ) -% -% PLOT_EXT.M 11.06.92 -% -% Plots the boundary and the extremal points: -% -% number = 1: extremal points are numbered. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% extract various variables from x - [length_t, t, alpha, r, lower_bound, lambda] = ... - def_var( param, x ); - - [t,I] = sort( t ); - - clg; - critical_points = param(20:length(param)); - t_old = t; - t = (0:.01:1)'; -% corners should be ploted in any case - for j=1:length(critical_points) - I = find( critical_points(j) == t ); - if length(I) == 0, t = [t; critical_points(j)]; end - end - t = sort(t); - - gamma = eval( BOUNDARYstr ); - -% compute visually pleasing axis coordinates - maRe = max( real( gamma ) ); - miRe = min( real( gamma ) ); - maIm = max( imag( gamma ) ); - miIm = min( imag( gamma ) ); - - - mRe = ( maRe + miRe )/2; - mIm = ( maIm + miIm )/2; - lRe = ( maRe - miRe ); - lIm = ( maIm - miIm ); - - if lRe > lIm - sg = 0.55 * lRe; - else - sg = 0.55 * lIm; - end %if - - -% plot boundary, extremal points - axis( [mRe-sg, mRe+sg, mIm-sg, mIm+sg] ); - axis( 'square' ); - plot(real(gamma),imag(gamma),'w'); - hold; - t = t_old; - gamma = eval( BOUNDARYstr ); - - plot( real(gamma),imag(gamma), 'ow' ); - plot([-100 100],[0 0],'w'); - plot([0 0],[-100 100],'w'); - - if number == 1 - for j=1:length_t - t = t_old(j); - gamma = eval( BOUNDARYstr ); - x = real( gamma ); - y = imag( gamma ); - text( x, y, int2str(j) ); - end % for - end; % if - - pause - hold - axis('normal') -return; //GO.SYSIN DD plot_ext.m echo q_newton.m 1>&2 sed >q_newton.m <<'//GO.SYSIN DD q_newton.m' 's/^-//' - function [ param, x, error_norm, relative_error ] = ... - q_newton( param, FUNstr, BASISstr, BOUNDARYstr, x ); -%function [ param, x, error_norm, relative_error ] = ... -%q_newton( param, FUNstr, BASISstr, BOUNDARYstr, x ); -% -% Q_NEWTON.M 05.07.92 -% -% Computes the solution to the nonlinear problem defined by -% NEWTON_F.M using a quasi Newton technique. -% For further information see SET_PARA.M. -% -% Copyright (C) 1992, Bernd Fischer, Jan Modersitzki. -% All rights reserved. - -% -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - R_dim = param(2); - stop_it = param(11); % if norm( search direction ) <= - % param(11) stop Newton iteration - delta = param(12); % stepsize for the numerical derivative - max_iter = param(14); % maximum number of newton iterations - output = param(17); - x_old = x; - - if output > 1 - disp( ['NEWTON'] ); - end; - - format long - if output > 0 - disp( [' #Iter norm( step ) norm( rhs ) ' ] ); - end; - - l = length(x); - length_t = (l - (R_dim + 1)) / 3; - G = zeros(l,l); - search_direction = ones(l,1); - iter = 1; - - - while norm( search_direction ) > stop_it & iter <= max_iter - - % compute approximation to the Hessian - for j=1:l - - old = x(j); - x(j) = x(j) + delta; - G(:,j) = newton_f( param, FUNstr, BASISstr, BOUNDARYstr, ... - x, length_t ); - - x(j) = old - delta; - G(:,j) = G(:,j) - ... - newton_f( param, FUNstr, BASISstr, BOUNDARYstr,... - x, length_t ); - x(j) = old; - - end % for - - % compute rigth hand side - rhs = - newton_f( param, FUNstr, BASISstr, BOUNDARYstr,... - x, length_t ) * 2 * delta; - - search_direction = G\rhs; - - if output > 0 - disp( [ iter, ... - norm( search_direction ) , ... - norm( rhs/( 2*delta) ) ] ); - end; - - x = x + search_direction; - iter = iter + 1; - - end %while - - format short - - if output > 3 - disp( ['End NEWTON' ] ); - end; - - param(13) = iter; - if norm( search_direction ) <= stop_it - param(15) = 9; % success - else - param(15) = 10; % no success - x = x_old; - end; - - % compute norm - [ error_norm, relative_error, t_max, alpha_max ] = ... - compnorm( param, FUNstr, BASISstr, BOUNDARYstr, x ); -return; - //GO.SYSIN DD q_newton.m echo recangle.m 1>&2 sed >recangle.m <<'//GO.SYSIN DD recangle.m' 's/^-//' - function [gamma, dgamma] = recangle( t, right, upper, left, lower ); -%function [gamma, dgamma] = recangle( t, right, upper, left, lower ); -% -% RECANGLE.M 16.07.92 -% -% The (real) parameters define a rectangle with corners -% (right + i*upper), (left +i*lower) in the complex plane. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - gamma = zeros( length(t), 1 ); - dgamma = ones( length(t), 1 ); - - I1 = find( t < 0.25 ); - if length(I1) > 0 - gamma(I1) = right + ... - i * ( lower + ( upper - lower ) * 4 * t(I1) ); - dgamma(I1) = i * ( upper - lower ) * 4 * dgamma(I1); - end; - - I2 = find( t >= 0.25 & t < 0.5 ); - if length(I2) > 0 - t2 = t(I2) - 0.25; - gamma(I2) = right - ( right - left ) * 4 * t2 + ... - i * upper; - dgamma(I2) = - ( right - left ) * 4 * dgamma(I2); - end; - - I3 = find( t >= 0.5 & t < 0.75 ); - if length(I3) > 0 - t3 = t(I3) - 0.5; - gamma(I3) = left +... - i * ( upper - ( upper - lower ) * 4 * t3 ); - dgamma(I3) = - i * ( upper - lower ) * 4 * dgamma(I3); - end; - - I = find(t >= 0.75 ); - if length(I) > 0 - t4 = t(I) - 0.75; - gamma(I) = left + ( right - left ) * 4 * t4 +... - i * lower; - dgamma(I) = ( right - left ) * 4 * dgamma(I); - end; -return //GO.SYSIN DD recangle.m echo remez.m 1>&2 sed >remez.m <<'//GO.SYSIN DD remez.m' 's/^-//' - function [ param, x, error_norm, relative_error ] = ... - remez( param, FUNstr, BASISstr, BOUNDARYstr, x ); -%function [ param, x, error_norm, relative_error ] = ... -%remez( param, FUNstr, BASISstr, BOUNDARYstr, x ); -% -% REMEZ.M 16.07.92 -% -% This function computes the coefficients lambda (part of x) of the -% best approximation (in terms of the basis functions defined in -% BASISstr) of the function (given by FUNstr) on a boundary -% (defined by BOUNDARYstr). As a by-product the function returns -% the extremal points (t, alpha), some weights r (related to the dual -% formulation of the problem), and a lower and upper bound for -% the minimal deviation. -% -% INPUT: -% param: vector of control parameters (see SET_PARA.M) -% FUNstr: name of approximant, (see CONVSTR.M) -% e.g.: FUNstr = 'fcos(gamma);' -% BASISstr: name of linear space, (see CONVSTR.M) -% e.g.: BASISstr = 'monom( gamma, C_dim, choice );' -% BOUNDARYstr: name of boundary function (see CONVSTR.M) -% e.g.: BOUNDARYstr = 'circle(t);' -% x: collection of variables: -% t actual extremal points -% alpha corresponding angles -% r weights -% lower_bound deviation at actual extremal points -% lambda actual coefficients of the approximation -% -% OUTPUT: -% param: actual parameters and communication switch -% x: [t; alpha; r; lower_bound; lambda] -% error_norm: upper bound for the minimal deviation -% relative_error: relative distance between upper and lower bound -% -% The function has three main ingredients: -% -% 1. INITIAL.M -% sets initial extremal points depending on param(3) -% (see: SET_EXTR.M), defines the initial system matrix A, -% and the initial right hand side c (see: SET_COLM.M). -% Moreover, the first set of weights r is computed. -% -% 2. BOUNDS.M -% computes upper and lower bounds for the minimal deviation -% and actual coefficients. -% Furthermore, it returns estimates for extremal points. -% -% 3. IMPROVE.M -% improves the actual approximation by exchanging 'old' -% extremal points with the 'new' ones obtained -% by BOUNDS.M. -% -% STOPITER.M: controls interrups -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - [ x, A, c ] = ... - initial( param, FUNstr, BASISstr, BOUNDARYstr, x ); - - [ x, error_norm, relative_error, t_max, alpha_max ] = ... - bounds( param, FUNstr, BASISstr, BOUNDARYstr, x, A, c ); - -% main loop - param(15) = stopiter( param, BOUNDARYstr, x, relative_error ); - - while param(15) == 0 % stopping criteria - - % iteration count - param(9) = param(9) + 1; - - if param(18) == 2 - picture( param, FUNstr, BASISstr, BOUNDARYstr, ... - x, error_norm, t_max, alpha_max ); - end; - - [ x, A, c ] = ... - improve( param, FUNstr, BASISstr, BOUNDARYstr,... - x, t_max, alpha_max, A, c ); - - [ x, error_norm, relative_error, t_max, alpha_max ] = ... - bounds( param, FUNstr, BASISstr, BOUNDARYstr, x, A, c ); - - param(15) = stopiter( param, BOUNDARYstr, x, relative_error ); - - end %while - - if param(18) == 2 - picture( param, FUNstr, BASISstr, BOUNDARYstr, ... - x, error_norm, t_max, alpha_max ); - end; - -% saveguard, saves actual data (see: SET_EXTR.M) - save remez -return; //GO.SYSIN DD remez.m echo sector.m 1>&2 sed >sector.m <<'//GO.SYSIN DD sector.m' 's/^-//' - function [gamma, dgamma] = sector( t, angle ); -%function [gamma, dgamma] = sector( t, angle ); -% -% SECTOR.M 16.07.92 -% -% t in [0,1], angle in [0, 360]. -% Defines an open circular sector in the complex plane. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - gamma = zeros( length(t), 1 ); - dgamma = ones( length(t), 1 ); - - - I1 = find( t <= 0.5 ); - if length(I1) > 0 - t1 = (angle * pi / 90) * t(I1); - gamma(I1) = cos(t1) + i*sin(t1); - dgamma(I1) = (angle * pi / 90) *... - ( - sin(t1) + i*cos(t1) ); - end; - - I = find(t > 0.5 ); - if length(I) > 0 - t2 = 2*(t(I) - 0.5); - gamma(I) = (1-t2) * exp( i * angle * 2 * pi / 360 ); - dgamma(I) = - 2 * exp( angle * 2 * pi / 90 )*... - ones(gamma(I)); - end; -return; //GO.SYSIN DD sector.m echo set_colm.m 1>&2 sed >set_colm.m <<'//GO.SYSIN DD set_colm.m' 's/^-//' - function [A,c] = set_colm( param, FUNstr, BASISstr, BOUNDARYstr, t, alpha ); -%function [A,c] = set_colm( param, FUNstr, BASISstr, BOUNDARYstr, t, alpha ); -% -% SET_COLM.M 16.07.92 -% -% Defines the m (= length(t)) columns of system matrix A -% and the right hand side c. -% -% Note: -% the j'th column of A, A(:,j) is defined by -% -% A(1,j) = 1; -% -% A(k,j) = real( exp( -i*alpha(j) ) * BASIS(k-1)(t(j))), -% k = 2,..., C_dim + 1, -% -% A(k,j) = real( exp( -i*alpha(j) ) * BASIS(k-C_dim-1)(t(j))), -% k = C_dim+2, ..., R_dim+1, -% -% where BASIS(k) denotes the k'th basis-function. -% -% For complex coefficients: -% basis-function BASIS(j+C_dim) := i*BASIS(j). -% -% c = real( exp( -i*alpha(j) ) * f(t(j)) ). -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - C_dim = param(1); - R_dim = param(2); - - [gamma] = eval( BOUNDARYstr ); - [phi] = eval( BASISstr ); - - if R_dim ~= C_dim - phi = [phi; i * phi]; - end - - m = length( t ); - A = ones( R_dim+1, m ); - A( 2:R_dim+1, : ) = real( phi*diag( exp( -i*alpha ) ) ); - - f = eval( FUNstr ); - c = real( exp(-i*alpha).*f ); -return; //GO.SYSIN DD set_colm.m echo set_extr.m 1>&2 sed >set_extr.m <<'//GO.SYSIN DD set_extr.m' 's/^-//' - function [t, alpha] = set_extr( param, x ); -%function [t, alpha] = set_extr( param, x ); -% -% SET_EXTR.M 16.07.92 -% -% Defines initial extremal points (t(i), alpha(i)) -% in [0,1]x[0,2pi], i = 1, ..., R_dim + 1. -% -% param(3) is used to control the initialization process: -% -% 1: equidistant -% 2: random -% 3: last run -% 4: input from a user supplied file -% 5: via keyboard -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - R_dim = param(2); - initial_extremal_points = param(3); - m = R_dim + 1; - param(3) = 0; % no further initialization - - if initial_extremal_points == 0 - - t = x(1:m); - alpha = x(m+1:2*m); - disp( ['old-values'] ); - - elseif initial_extremal_points == 1 % equidistant - - t = 0:1/R_dim:1; - alpha = pi * ( 2 * t - 1 ); - - elseif initial_extremal_points == 2 % random - - t = rand( m, 1 ); - alpha = rand( m, 1 ); - alpha = -pi + 2 * pi * alpha; - - elseif initial_extremal_points == 3 % last run - - load remez - t = x(1:m); - alpha = x(m+1:2*m); - - elseif initial_extremal_points == 4 % via file - - file=input('file: ','s'); - eval(['load ',file]) - - elseif initial_extremal_points == 5 % via keyboard - - clc; - t=[]; - alpha=[]; - while length(t) ~= m - disp( ['enter ',num2str(m),' numbers t_j in [0,1]'] ); - t = input( '[t_1 t_2 ... t_m] : ' ); - end; - - while length(alpha) ~= m - disp(['enter ',num2str(m),' numbers alpha_j in [-pi,pi]']); - alpha=input('[alpha_1 alpha_2 ... alpha_m] : '); - end; - end; %if - -% final shape - [dummy,col] = size(t); if col > 1, t=t'; end; - [dummy,col] = size(alpha); if col > 1, alpha=alpha'; end; - t = sort( t ); -return; //GO.SYSIN DD set_extr.m echo set_para.m 1>&2 sed >set_para.m <<'//GO.SYSIN DD set_para.m' 's/^-//' - function [param] = set_para( critical_points, C_dim, real_coef ); -%function [param] = set_para( critical_points, C_dim, real_coef ); -% -% SET_PARA.M 16.07.92 -% -% This routine fills the parameter vector param with default -% values and resets the variable x. -% -% INPUT: -% critical_points: parameter values of critical points on the -% boundary, e.g., corners of polygons etc., -% see also Q_NEWTON.M -% C_dim: C-dimension of the approximation space -% real_coef: if real_coef = 1, the coefficients are real -% -% OUTPUT: -% param(1) C-Dimension -% param(2) R-Dimension -% param(3) definition of initial extremal points (see: SET_EXTR.M) -% param(4) 1: single exchange step in IMPROVE.M -% 0: multiple exchange in IMPROVE.M -% param(5) stepsize in COMPNORM.M is 1/param(5) -% param(6) bound for relative error (see: STOPITER.M) -% param(7) if the relative error is below param(7), -% the algorithm checks for clustered extremal points -% (see: STOPITER.M) -% param(8) if the distance between extremal points is less then -% param(8), the user is prompted for clustering -% (see: STOPITER.M) -% param(9) counts iteration of remez (see REMEZ.M) -% param(10) maximum number of iterations -% param(11) if norm( search direction ) <= param(11) stop Newton iteration -% (see: Q_NEWTON.M) -% param(12) stepsize for the numerical derivative (see: Q_NEWTON.M) -% param(13) counts newton iterations (see: Q_NEWTON.M) -% param(14) maximum number of newton iterations -% param(15) return value, controls interrups (see: STOPMENU.M) -% param(16) used by the algorithm for controling events (see: MAIN.M) -% -% param(17) controls amount of output: -% 0: no output during iteration -% 1: number of iterations, error norm, lower bound, -% relative error -% 2: 1 + name of executed function -% 3: 2 + C_dim, R_dim, t, alpha, r, c -% t_max, alpha_max, error_max -% 4: 3 + pause -% -% param(18) plots: -% 1: after last remez iteration -% 2: after each remez iteration -% param(21,length(param)): -% critical points on the boundary -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. - -% -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - param = zeros( 20,1 ); - - param(1) = C_dim; - - if real_coef - param(2) = C_dim; - else - param(2) = 2 * C_dim; - end; - - param(3) = 2; - param(4) = 1; - param(5) = 50; - param(6) = 1e-3; - param(7) = 1e-5; - param(8) = 1E-4; - param(9) = 0; - param(10) = 100; - param(11) = 1.E-10; - param(12) = 1.E-5; - param(13) = 0; - param(14) = 10; - param(15) = 0; - param(16) = 1; - param(17) = 1; - param(18) = 0; - - param = [param; critical_points]; -return; //GO.SYSIN DD set_para.m echo sign_spe.m 1>&2 sed >sign_spe.m <<'//GO.SYSIN DD sign_spe.m' 's/^-//' - function [sign_spez] = sign_spe( z ); -%function [sign_spez] = sign_spe( z ); -% -% SIGN_SPE.M 16.07.92 -% -% Special signum function for usage in COMPNORM.M -% -% sign_spe := -1, x <= 0, -% 1, x > 0. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - sign_spez = ones(z); - I = find( z <= 0 ); - sign_spez(I) = -1 * ones(I); -return; //GO.SYSIN DD sign_spe.m echo simplex.m 1>&2 sed >simplex.m <<'//GO.SYSIN DD simplex.m' 's/^-//' - function [x, I] = simplex( c, A, b, I ); -%function [x, I] = simplex( c, A, b, I ); -% -% SIMPLEX.M 16.07.92 -% -% Simplex algorithm for solving -% -% c'x = min! -% subject to A*x = b, -% x >= 0, -% -% with c, x in R^n, b in R^m, and A in R^(m,n). -% -% Note: the algorithm assumes the knowledge of a start basis, defined -% by the index set I. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - -% No part of this code may be reproduced, stored in a retrieval system, -% translated, transcribed, transmitted, or distributed in any form -% or by any means, means manual, electric, electronic, electro-magnetic, -% mechanical, chemical, optical, photocopying, recording, or otherwise, -% without the prior explicit written permission of the authors or their -% designated proxies. In no event shall the above copyright notice be -% removed or altered in any way. -% -% This code is provided "as is", without any warranty of any kind, either -% expressed or implied, including but not limited to, any implied warranty -% of merchantibility or fitness for any purpose. In no event will any party -% who distributed the code be liable for damages or for any claim(s) by -% any other party, including but not limited to, any lost profits, lost -% monies, lost data or data rendered inaccurate, losses sustained by -% third parties, or any other special, incidental or consequential damages -% arrising out of the use or inability to use the program, even if the -% possibility of such damages has been advised against. The entire risk -% as to the quality, the performace, and the fitness of the program for any -% particular purpose lies with the party using the code. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Any use of this code constitutes acceptance of the terms of the above -% statements -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - x = zeros(b); - - end_loop = 0; - while end_loop == 0 - - % complement of I - J = []; - for j = 1:length(c) - if length( find( j == I ) ) == 0 - J = [J,j]; - end; % if - end; % for - - % find exchange candidate in J - B = A(:,I); - y = B'\c(I); - s = A'* y; - t = c(J) - s(J); - [t_r,r] = min(t); - r = J(r); - - if t_r >= -1e-15 % x is solution - % Note: the obvious request t_r >= 0 - % may cause trouble due to roundoff errors - - end_loop = 1; - - else - - % find exchange candidate in I - d_r = B\A(:,r); - x = B\b; - K = find( d_r > 0 ); - - if length(K) == 0 % objective function unbounded - - end_loop = 2; - - else - - % exchange - dummy = x(K)./ d_r(K); - [delta,q] = min( dummy ); - q = I(K(q)); - I = [I(find( q ~= I )), r]; - - end; %if - end; %if - end; %while - - if end_loop == 2 - - error( ['objective function unbounded!'] ); - - else - - % compute solution - z = zeros(c); - x = B\b; - z(I) = x; - x = z; - - end; %if -return; - - //GO.SYSIN DD simplex.m echo stopiter.m 1>&2 sed >stopiter.m <<'//GO.SYSIN DD stopiter.m' 's/^-//' - function [value] = ... - stopiter( param, BOUNDARYstr, x, relative_error ); -%function [value] = ... -%stopiter( param, BOUNDARYstr, x, relative_error ); -% -% STOPITER.M 16.07.92 -% -% Dependent on the input parameter and the program status, this -% function returns an interrupt value: -% -% 0: continue, -% 1: relative_error <= error_bound, -% 2: number_iter > max_iter, -% 4: dist(extremal points) <= coincide_bound. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - - m = param(2) + 1; - error_bound = param(6); - cluster_bound = param(7); - coincide_bound = param(8); - number_iter = param(9); - max_iter = param(10); - newton_iter = param(14); - t = x(1:m); - - if relative_error <= error_bound - err = 1; - else - err = 0; - end; - - if number_iter >= max_iter - iter = 2; - else - iter = 0; - end; - - if relative_error <= cluster_bound & newton_iter > 0 - - t = sort(t); - - if min( abs( t(2:length(t) ) - t(1:length(t)-1) ) ) > coincide_bound - - ext = 0; - - else - - number = 1; % number extremal points in PLOT_EXT.M - plot_ext( param, BOUNDARYstr, x, number ) - disp(['relative error = ', num2str(relative_error) ]); - str = input('Do you want to cluster [n] ', 's'); - - if str == 'y' - ext = 4; - else - ext = 0; - end; % if - - end; % if - - else - - ext = 0; - - end; % if - - value = ext + iter + err; -return //GO.SYSIN DD stopiter.m echo stopmenu.m 1>&2 sed >stopmenu.m <<'//GO.SYSIN DD stopmenu.m' 's/^-//' - function [param] = stopmenu( param, rel_error ); -%function [param] = stopmenu( param, rel_error ); -% -% STOPMENU.M 16.07.92 -% -% param(15): input state of program -% param(16): output event loop (see COCA.M) -% 1 : call remez -% 2 : call newton -% 3 : cluster extremal-points -% 4 : show extremal-points -% 5 : show error-curves -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - clc - - m0 = 'Exit PROGRAM'; - m1 = 'Call REMEZ with new error bound'; - m2 = 'Call REMEZ with new max of iterations'; - m3 = 'Call REMEZ with new max of iter & err bound'; - m4 = 'Cluster extremal-points'; - m5 = 'Call NEWTON'; - m6 = 'Show extremal-points'; - m7 = 'Show error-curves'; - m8 = 'Call REMEZ'; - - if param(15) == 1 % relative error <= param(6) - disp( ['relative error = ', num2str( rel_error ), ... - ' <= ', num2str(param(6)) ] ); - - if param(14) ~= 0 - Mstr = menustr( 'REMEZ: success', ... - m0, m1, m5, m6, m7 ); - else - Mstr = menustr( 'REMEZ: success', ... - m0, m1, m6, m7 ); - end - - elseif param(15) == 2 % #iter > max_iter - disp( ['number of iteration > ', int2str(param(10)) ] ); - - if param(14) ~= 0 - Mstr = menustr( 'REMEZ: iterations', m0, m2, m5, m6, m7 ) ; - else - Mstr = menustr( 'REMEZ: iterations', m0, m2, m6, m7 ) ; - end; - - elseif param(15) == 3 % #iter > max_iter - % rel error <= error bound - disp( ['relative error <= ', num2str(param(6)) ] ); - disp( ['number of iteration > ', int2str(param(10)) ] ); - - if param(14) ~= 0 - Mstr = menustr( 'REMEZ: success, iteration', ... - m0, m3, m5, m6, m7 ); - else - Mstr = menustr( 'REMEZ: success, iteration', ... - m0, m3, m6, m7 ); - end; - - elseif param(15) == 4 % extremal-points coincide - Mstr= m4; - - elseif param(15) == 5 % extremal-points coincide - % rel error <= error bound - disp( ['relative error <= ', num2str(param(6)) ] ); - Mstr = menustr( 'REMEZ: success, cluster', ... - m0, m1, m4, m6, m7) ; - - elseif param(15) == 6 % #iter > max_iter - % extremal-points coincide - disp( ['number of iteration > ', int2str(param(10)) ] ); - Mstr = menustr( 'REMEZ: iterations, cluster', ... - m0, m2, m4, m6, m7) ; - - elseif param(15) == 7 % #iter > max_iter - % extremal-points coincide - % rel error <= error bound - disp( ['relative error <= ', num2str(param(6)) ] ); - disp( ['number of iteration > ', int2str(param(10)) ] ); - Mstr = menustr( 'REMEZ: success, iterations, cluster',... - m0, m3, m4, m6, m7) ; - - elseif param(15) == 8 % back from clustering - Mstr=m5; - - elseif param(15) == 9 % newton successful terminated - Mstr = menustr( 'Newton: success', m0, m6, m7) ; - - elseif param(15) == 10 % newton without success - Mstr = menustr( 'Newton: no success', m0, m8, m6, m7) ; - - end; - - if strcmp( Mstr, m0 ) - param(16) = 0; - - elseif strcmp( Mstr, m1 ) - param(3) = 0; - param(6) = input( 'new error bound: ' ); - param(15) = 0; - param(16) = 1; - - elseif strcmp( Mstr, m2 ) - param(3) = 0; - param(10) = input( 'new max iterations: ' ); - param(15) = 0; - param(16) = 1; - - elseif strcmp( Mstr, m3 ) - param(3) = 0; - param(6) = input( 'new error bound: ' ); - param(10) = input( 'new max iterations: ' ); - param(15) = 0; - param(16) = 1; - - elseif strcmp( Mstr, m4 ) - param(16) = 3; - - elseif strcmp( Mstr, m5 ) - param(16) = 2; - - elseif strcmp( Mstr, m6 ) - param(16) = 4; - - elseif strcmp( Mstr, m7 ) - param(16) = 5; - - elseif strcmp( Mstr, m8 ) - param(15) = 0; - param(16) = 1; - - else - error('---') - end; -return; //GO.SYSIN DD stopmenu.m echo two_circ.m 1>&2 sed >two_circ.m <<'//GO.SYSIN DD two_circ.m' 's/^-//' - function [gamma, dgamma] = two_circ( t, Rem1, Imm1, r1, Rem2, Imm2, r2 ); -%function [gamma, dgamma] = two_circ( t, Rem1, Imm1, r1, Rem2, Imm2, r2 ); -% -% TWO_CIRC.M 16.07.92 -% -% Defines a boundary consisting of two circles, one with center -% (Rem1 + i*Imm1) and radius r1, the other with center -% (Rem2 + i*Imm2) and radius r2. -% -% This function is part of the COCA - package. -% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki. -% All rights reserved. -% - m1 = Rem1 + i * Imm1; - m2 = Rem2 + i * Imm2; - - gamma = zeros(t); - dgamma = zeros(t); - - I1 = find( t < 0.5 ); - if length(I1) > 0 - t1 = 4*pi* t(I1); - gamma(I1) = r1 * (cos(t1) + i*sin(t1)) + m1; - dgamma(I1) = (-sin(t1) + i*cos(t1))*4*pi*r1; - end; - - I2 = find( t >= 0.5 ); - if length(I2) > 0 - t2 = 4*pi*(t(I2) - 0.5); - gamma(I2) = r2 * (cos(t2) + i*sin(t2)) + m2; - dgamma(I2) = (-sin(t2) + i*cos(t2) )*4*pi*r2; - end; -return //GO.SYSIN DD two_circ.m .