1 ## Copyright (C) 2011 Fernando Damian Nieuwveldt <fdnieuwveldt@gmail.com>
2 ## 2012 Adapted by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
4 ## This program is free software; you can redistribute it and/or
5 ## modify it under the terms of the GNU General Public License
6 ## as published by the Free Software Foundation; either version 3
7 ## of the License, or (at your option) any later version.
9 ## This program is distributed in the hope that it will be useful,
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
14 ## @deftypefn {Function File} {} cauchy (@var{N}, @var{r}, @var{x}, @var{f} )
15 ## Return the Taylor coefficients and numerical differentiation of a function
16 ## @var{f} for the first @var{N-1} coefficients or derivatives using the fft.
17 ## @var{N} is the number of points to evaluate,
18 ## @var{r} is the radius of convergence, needs to be chosen less then the smallest singularity,
19 ## @var{x} is point to evaluate the Taylor expansion or differentiation. For example,
21 ## If @var{x} is a scalar, the function @var{f} is evaluated in a row vector
22 ## of length @var{N}. If @var{x} is a column vector, @var{f} is evaluated in a
23 ## matrix of length(x)-by-N elements and must return a matrix of the same size.
27 ## d = cauchy(16, 1.5, 0, @@(x) exp(x));
28 ## @result{} d(2) = 1.0000 # first (2-1) derivative of function f (index starts from zero)
33 function deriv = cauchy(N, r, x, f)
41 error('cauchy:InvalidArgument', 'The 3rd argument must be a column vector');
47 f_p = f (bsxfun (@plus, x, r * exp (i * th) ) );
49 evalfft = real(fft (f_p, [], 2));
51 deriv = bsxfun (@times, evalfft, 1./(N*(r.^n)).* factorial(n)) ;
55 function g = hermite(order,x)
56 ## N should be bigger than order+1
59 Hnx = @(t) exp ( bsxfun (@minus, kron(t(:).', x(:)) , t(:).'.^2/2) );
60 Hnxfft = cauchy(N, r, 0, Hnx);
61 g = Hnxfft(:, order+1);
65 %! # Cauchy integral formula: Application to Hermite polynomials
66 %! # Author: Fernando Damian Nieuwveldt
67 %! # Edited by: Juan Pablo Carbajal
69 %! Hnx = @(t,x) exp ( bsxfun (@minus, kron(t(:).', x(:)) , t(:).'.^2/2) );
70 %! hermite = @(order,x) cauchy(32, 0.5, 0, @(t)Hnx(t,x))(:, order+1);
72 %! t = linspace(-1,1,30);
73 %! he2 = hermite(2,t);
78 %! plot(t,he2,'bo;Contour integral representation;', t,he2_,'r;Exact;');
82 %! % --------------------------------------------------------------------------
83 %! % The plots compares the approximation of the Hermite polynomial using the
84 %! % Cauchy integral (circles) and the corresposind polynomial H_2(x) = x.^2 - 1.
85 %! % See http://en.wikipedia.org/wiki/Hermite_polynomials#Contour_integral_representation
88 %! # Cauchy integral formula: Application to Hermite polynomials
89 %! # Author: Fernando Damian Nieuwveldt
90 %! # Edited by: Juan Pablo Carbajal
92 %! xx = sort (rand (100,1));
93 %! yy = sin (3*2*pi*xx);
95 %! # Exact first derivative derivative
96 %! diffy = 6*pi*cos (3*2*pi*xx);
98 %! np = [10 15 30 100];
101 %! idx = sort(randperm (100,np(i)));
106 %! yval = ppval (ppder(p),x);
107 %! # Use the cauchy formula for computing the derivatives
108 %! deriv = cauchy (fix (np(i)/4), .1, x, @(x) sin (3*2*pi*x));
111 %! h = plot(xx,diffy,'-b;Exact;',...
112 %! x,yval,'-or;ppder solution;',...
113 %! x,deriv(:,2),'-og;Cauchy formula;');
114 %! set(h(1),'linewidth',2);
115 %! set(h(2:3),'markersize',3);
117 %! legend(h, 'Location','Northoutside','Orientation','horizontal');
123 %! % --------------------------------------------------------------------------
124 %! % The plots compares the derivatives calculated with Cauchy and with ppder.
125 %! % Each subplot shows the results with increasing number of samples.