function c = fdcoeffF(k,xbar,x) % Compute coefficients for finite difference approximation for the % derivative of order k at xbar based on grid values at points in x. % % This function returns a row vector c of dimension 1 by n, where n=length(x), % containing coefficients to approximate u^{(k)}(xbar), % the k'th derivative of u evaluated at xbar, based on n values % of u at x(1), x(2), ... x(n). % % If U is a column vector containing u(x) at these n points, then % c*U will give the approximation to u^{(k)}(xbar). % % Note for k=0 this can be used to evaluate the interpolating polynomial % itself. % % Requires length(x) > k. % Usually the elements x(i) are monotonically increasing % and x(1) <= xbar <= x(n), but neither condition is required. % The x values need not be equally spaced but must be distinct. % % This program should give the same results as fdcoeffV.m, but for large % values of n is much more stable numerically. % % Based on the program "weights" in % B. Fornberg, "Calculation of weights in finite difference formulas", % SIAM Review 40 (1998), pp. 685-691. % % Note: Forberg's algorithm can be used to simultaneously compute the % coefficients for derivatives of order 0, 1, ..., m where m <= n-1. % This gives a coefficient matrix C(1:n,1:m) whose k'th column gives % the coefficients for the k'th derivative. % % In this version we set m=k and only compute the coefficients for % derivatives of order up to order k, and then return only the k'th column % of the resulting C matrix (converted to a row vector). % This routine is then compatible with fdcoeffV. % It can be easily modified to return the whole array if desired. % % From http://www.amath.washington.edu/~rjl/fdmbook/ (2007) n = length(x); if k >= n error('*** length(x) must be larger than k') end m = k; % change to m=n-1 if you want to compute coefficients for all % possible derivatives. Then modify to output all of C. c1 = 1; c4 = x(1) - xbar; C = zeros(n-1,m+1); C(1,1) = 1; for i=1:n-1 i1 = i+1; mn = min(i,m); c2 = 1; c5 = c4; c4 = x(i1) - xbar; for j=0:i-1 j1 = j+1; c3 = x(i1) - x(j1); c2 = c2*c3; if j==i-1 for s=mn:-1:1 s1 = s+1; C(i1,s1) = c1*(s*C(i1-1,s1-1) - c5*C(i1-1,s1))/c2; end C(i1,1) = -c1*c5*C(i1-1,1)/c2; end for s=mn:-1:1 s1 = s+1; C(j1,s1) = (c4*C(j1,s1) - s*C(j1,s1-1))/c3; end C(j1,1) = c4*C(j1,1)/c3; end c1 = c2; end c = C(:,end)'; % last column of c gives desired row vector