-
Notifications
You must be signed in to change notification settings - Fork 29
/
numdiff.m
78 lines (61 loc) · 3.16 KB
/
numdiff.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
function [df, d2f]=numdiff(x,f)
% NUMDIFF : compute first and second order derivative using possibly unequally spaced
% data.
%
% Syntax
% [df d2f] = numdiff(x, f)
%
% Description
% numdiff computes first and second order derivatives using a 5th order
% formula that accounts for possibly unequally spaced data. Because a
% 5th order centered difference formula is used, numdiff returns NaNs for
% the first 2 and last 2 points in the input vector for x.
%
% Input :
% x - a vector of values of x
% f - a vector of values of the function f for which derivatives are to
% be computed. Must be the same length as x.
%
% Output:
% df - a vector of length(x) containing the first derivative of f at
% each point x except at the first 2 and last 2 points.
% d2f - a vector of length(x) containing the second derivative of f at
% each point x except at the first 2 and last 2 points.
%
% Sources:
% [1] M. K. Bowen, R. Smith, "Derivative formulae and errors for
% non-uniformly spaced points", Proceedings of the Royal Society A,
% vol. 461 pp 1975-1997, July 2005. DOI: 10.1098/rpsa.2004.1430
% convert to column vectors if needed
f=f(:);
x=x(:);
N=length(f);
df = NaN(size(f));
% first two points are special
df(1:2) = NaN;
% last two points are special
df((N-1):N) = NaN;
% rest of points. Take reference point to be the middle of each group of 5
% points. Calculate displacements
F = [f(1:(N-4)) f(2:(N-3)) f(3:(N-2)) f(4:(N-1)) f(5:N)];
A = [x(1:(N-4)) x(2:(N-3)) x(3:(N-2)) x(4:(N-1)) x(5:N)] - repmat(x(3:(N-2)),1,5);
% A(:,3) is always zero but we retain it in case we want to generalize this
% function in the future
U(:,1) = A(:,2).*A(:,3).*A(:,4) + A(:,2).*A(:,3).*A(:,5) + A(:,2).*A(:,4).*A(:,5) + A(:,3).*A(:,4).*A(:,5);
U(:,2) = A(:,1).*A(:,3).*A(:,4) + A(:,1).*A(:,3).*A(:,5) + A(:,1).*A(:,4).*A(:,5) + A(:,3).*A(:,4).*A(:,5);
U(:,3) = A(:,1).*A(:,2).*A(:,4) + A(:,1).*A(:,2).*A(:,5) + A(:,1).*A(:,4).*A(:,5) + A(:,2).*A(:,4).*A(:,5);
U(:,4) = A(:,1).*A(:,2).*A(:,3) + A(:,1).*A(:,2).*A(:,5) + A(:,1).*A(:,3).*A(:,5) + A(:,2).*A(:,3).*A(:,5);
U(:,5) = A(:,1).*A(:,2).*A(:,3) + A(:,1).*A(:,2).*A(:,4) + A(:,1).*A(:,3).*A(:,4) + A(:,2).*A(:,3).*A(:,4);
L(:,1) = (A(:,1) - A(:,2)).*(A(:,1) - A(:,3)).*(A(:,1) - A(:,4)).*(A(:,1) - A(:,5));
L(:,2) = (A(:,2) - A(:,1)).*(A(:,2) - A(:,3)).*(A(:,2) - A(:,4)).*(A(:,2) - A(:,5));
L(:,3) = (A(:,3) - A(:,1)).*(A(:,3) - A(:,2)).*(A(:,3) - A(:,4)).*(A(:,3) - A(:,5));
L(:,4) = (A(:,4) - A(:,1)).*(A(:,4) - A(:,2)).*(A(:,4) - A(:,3)).*(A(:,4) - A(:,5));
L(:,5) = (A(:,5) - A(:,1)).*(A(:,5) - A(:,2)).*(A(:,5) - A(:,3)).*(A(:,5) - A(:,4));
df(3:(N-2)) = sum(-(U./L).*F,2);
% second derivative
U2(:,1) = A(:,2).*A(:,3)+A(:,2).*A(:,4)+A(:,2).*A(:,5)+A(:,3).*A(:,4)+A(:,3).*A(:,5)+A(:,4).*A(:,5);
U2(:,2) = A(:,1).*A(:,3)+A(:,1).*A(:,4)+A(:,1).*A(:,5)+A(:,3).*A(:,4)+A(:,3).*A(:,5)+A(:,4).*A(:,5);
U2(:,3) = A(:,1).*A(:,2)+A(:,1).*A(:,4)+A(:,1).*A(:,5)+A(:,2).*A(:,4)+A(:,2).*A(:,4)+A(:,4).*A(:,5);
U2(:,4) = A(:,1).*A(:,2)+A(:,1).*A(:,3)+A(:,1).*A(:,5)+A(:,2).*A(:,3)+A(:,2).*A(:,5)+A(:,3).*A(:,5);
U2(:,5) = A(:,1).*A(:,2)+A(:,1).*A(:,3)+A(:,1).*A(:,4)+A(:,2).*A(:,3)+A(:,2).*A(:,5)+A(:,3).*A(:,4);
d2f(3:(N-2)) = 2*sum(U2.*F,2);