function y = FSFun(x) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % y = FSFun(x) % ____________________________________________________________________________ % 18.385 (Nonlinear Dynamical Systems). MIT, Fall 1999. R. R. Rosales. % % User provided function for input into fourierSC or heatSln. Arbitrary % periodic functions F = F(x) --- of period 2*pi --- are allowed, provided a % file FSFun.m (with a function script) is present in the directory from which % the user is running MatLab. The MatLab statement y = FSFun(x) should return % the values of F(x) in the array y [i.e. y(n)=F(x(n))], for any array x with % x(n) in [0, 2*pi]. The function in FSFun will be used if the appropriate % option picked when running fourierSC or heatSln (a default is provided). % ^^^^^^^^^^^^^^^^^^^^^ % A simple example of a FSFun.m file is: % % function y = FSFun(x) % y = exp(-10*(x-pi).^2); % % Note the use of the vector operand .^, needed to get a vector valued answer. % % NOTE: the function whose Fourier series is being computed is ASSUMED to be % piecewise continuous. The discontinuities must be separated by a % distance of at least 2*pi/50 with the current resolution in the scripts. The % number of them cannot be too large either. The evaluation near a point of % discontinuity --- say, near x_d --- must be done in such a way that % % (A) F(x) defined by continuous expressions for both x > x_d and x < x_d; % with some value assigned at x = x_d. % (B) The value assigned at x_d satisfies: % % F(x_d) = Fl or F(x_d) = Fr or F(x_d) = 0.5*(Fl + Fr), % % where: Fl and Fr are the limits of F(x) from the left and right at x_d. % These three alternatives are useful in preserving nice symmetry (odd or % even) properties of the function F when evaluating it on a discrete set % of points x(j). % % For example, a "square wave" may be evaluated by the script % % function y = FSFun(x) % y = ones(1, length(x)); % index = find((x > 0.5) & (x < 2.5)); % y(index) = zeros(1, length(index)); % % 1) The first statement here creates an array of ones of the same size as x. % 2) The second finds the sets of indexes j such that 0.5 < x(j) < 2.5. % 3) The third statement sets the value of y to zero at the j's fund in (2). % % A special value could be assigned at the discontinuity points. For example, % add the statements % % Idis = find(x == 0.5); % y(Idis) = CONSTANT; % % where CONSTANT = 0 or 0.5 or 1 (here F+ = 1 and F- = 0). Note that Idis may % be the empty array if there is no j for which x(j) = 0.5. Then the second % statement here does not affect y. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fprintf('\n NOTE: the default choice for FSFun is being used.') fprintf('\n Do help FSFun to see how to use a different function.\n') % % y = (x.^2).*(x-2*pi).^2; % EXAMPLE with discontinuous 2nd derivatives. % % ------ Example with asymmetric square wave with discontinuity close to edge % of domain. Used to test part of code that creates the refined grid % xRef in fourierSC.m y = zeros(size(x)); n1 = min(10, ceil(length(x)/10)); n2 = ceil(length(x)/4); y(n1:(n1+n2)) = ones(1, n2+1); % % EOF