# to unbundle, sh this file (in an empty directory)
echo README 1>&2
-
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
-% 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
-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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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 <y|n> [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
- 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
- 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
.
