File Coverage

blib/lib/Math/Derivative.pm
Criterion Covered Total %
statement 56 59 94.9
branch 5 10 50.0
condition n/a
subroutine 8 8 100.0
pod 2 3 66.6
total 71 80 88.7


line stmt bran cond sub pod time code
1             package Math::Derivative;
2              
3 6     6   51722 use 5.010001;
  6         22  
4 6     6   31 use Exporter;
  6         9  
  6         522  
5              
6             our @ISA = qw(Exporter);
7             our %EXPORT_TAGS = (all => [qw(
8             Derivative1
9             Derivative2
10             centraldiff
11             forwarddiff
12             )]);
13              
14             our @EXPORT_OK = (@{$EXPORT_TAGS{all}});
15              
16             our $VERSION = 1.01;
17              
18 6     6   35 use strict;
  6         19  
  6         185  
19 6     6   35 use warnings;
  6         17  
  6         163  
20 6     6   31 use Carp;
  6         13  
  6         3897  
21              
22             =head1 NAME
23              
24             Math::Derivative - Numeric 1st and 2nd order differentiation
25              
26             =head1 SYNOPSIS
27              
28             use Math::Derivative qw(:all);
29              
30             @dydx = forwarddiff(\@x, \@y);
31              
32             @dydx = centraldiff(\@x, \@y);
33              
34             @dydx = Derivative1(\@x, \@y); # A synonym for centraldiff()
35              
36             @d2ydx2 = Derivative2(\@x, \@y, $yd0, $ydn);
37              
38             =head1 DESCRIPTION
39              
40             This Perl package exports functions that numerically approximate first
41             and second order differentiation on vectors of data. The accuracy of
42             the approximation will depend upon the differences between the
43             successive values in the X array.
44              
45             =head2 FUNCTIONS
46              
47             The functions may be imported by name or by using the tag ":all".
48              
49             =head3 forwarddiff()
50              
51             @dydx = forwarddiff(\@x, \@y);
52              
53             Take the references to two arrays containing the x and y ordinates of
54             the data, and return an array of approximate first derivatives at the
55             given x ordinates, using the forward difference approximation.
56              
57             The last term is actually formed using a backward difference formula,
58             there being no array item to subtract from at the end of the array.
59             If you want to use derivatives strictly formed from the forward
60             difference formula, use only the values from [0 .. #y-1], e.g.:
61              
62             @dydx = (forwarddiff(\@x, \@y))[0 .. $#y-1];
63              
64             or, more simply,
65              
66             @dydx = forwarddiff(\@x, \@y);
67             pop @dydx;
68              
69             =cut
70              
71             sub forwarddiff
72             {
73 1     1 1 80 my($x, $y) = @_;
74 1         1 my @y2;
75 1         2 my $n = $#{$x};
  1         2  
76              
77 1 50       1 croak "X and Y array lengths don't match." unless ($n == $#{$y});
  1         3  
78              
79 1         6 $y2[$n] = ($y->[$n] - $y->[$n-1])/($x->[$n] - $x->[$n-1]);
80              
81 1         4 for my $i (0 .. $n-1)
82             {
83 8         13 $y2[$i] = ($y->[$i+1] - $y->[$i])/($x->[$i+1] - $x->[$i]);
84             }
85              
86 1         4 return @y2;
87             }
88              
89             =head3 centraldiff()
90              
91             @dydx = centraldiff(\@x, \@y);
92              
93             Take the references to two arrays containing the x and y ordinates of
94             the data, and return an array of approximate first derivatives at the
95             given x ordinates.
96              
97             The algorithm used three data points to calculate the derivative, except
98             at the end points, where by necessity the forward difference algorithm
99             is used instead. If you want to use derivatives strictly formed from
100             the central difference formula, use only the values from [1 .. #y-1],
101             e.g.:
102              
103             @dydx = (centraldiff(\@x, \@y))[1 .. $#y-1];
104              
105             =cut
106              
107             sub centraldiff
108             {
109 3     3 1 256 my($x, $y) = @_;
110 3         5 my @y2;
111 3         6 my $n = $#{$x};
  3         6  
112              
113 3 50       6 croak "X and Y array lengths don't match." unless ($n == $#{$y});
  3         8  
114              
115 3         15 $y2[0] = ($y->[1] - $y->[0])/($x->[1] - $x->[0]);
116 3         14 $y2[$n] = ($y->[$n] - $y->[$n-1])/($x->[$n] - $x->[$n-1]);
117              
118 3         11 for my $i (1 .. $n-1)
119             {
120 12         29 $y2[$i] = ($y->[$i+1] - $y->[$i-1])/($x->[$i+1] - $x->[$i-1]);
121             }
122              
123 3         11 return @y2;
124             }
125              
126             =head3 Derivative2()
127              
128             @d2ydx2 = Derivative2(\@x, \@y);
129              
130             or
131              
132             @d2ydx2 = Derivative2(\@x, \@y, $yp0, $ypn);
133              
134             Take references to two arrays containing the x and y ordinates of the
135             data and return an array of approximate second derivatives at the given
136             x ordinates.
137              
138             You may optionally give values to use as the first derivatives at the
139             start and end points of the data. If you don't, first derivative values
140             will be assumed to be zero.
141              
142             =cut
143              
144             sub seconddx
145             {
146 1     1 0 354 my($x, $y, $yp1, $ypn) = @_;
147 1         3 my(@y2, @u);
148 1         1 my $n = $#{$x};
  1         2  
149              
150 1 50       2 croak "X and Y array lengths don't match." unless ($n == $#{$y});
  1         3  
151              
152 1 50       3 if (defined $yp1)
153             {
154 1         2 $y2[0] = -0.5;
155 1         4 $u[0] = (3/($x->[1] - $x->[0])) *
156             (($y->[1] - $y->[0])/($x->[1] - $x->[0]) - $yp1);
157             }
158             else
159             {
160 0         0 $y2[0] = 0;
161 0         0 $u[0] = 0;
162             }
163              
164 1         3 for my $i (1 .. $n-1)
165             {
166 23         35 my $sig = ($x->[$i] - $x->[$i-1])/($x->[$i+1] - $x->[$i-1]);
167 23         32 my $p = $sig * $y2[$i-1] + 2.0;
168              
169 23         29 $y2[$i] = ($sig - 1.0)/$p;
170 23         56 $u[$i] = (6.0 * (
171             ($y->[$i+1] - $y->[$i])/($x->[$i+1] - $x->[$i]) -
172             ($y->[$i] - $y->[$i-1])/($x->[$i] - $x->[$i-1]))/
173             ($x->[$i+1] - $x->[$i-1]) - $sig * $u[$i-1])/$p;
174             }
175              
176 1 50       3 if (defined $ypn)
177             {
178 1         1 my $qn = 0.5;
179 1         3 my $un = (3.0/($x->[$n]-$x->[$n-1])) *
180             ($ypn - ($y->[$n] - $y->[$n-1])/($x->[$n] - $x->[$n-1]));
181 1         2 $y2[$n] = ($un - $qn * $u[$n-1])/($qn * $y2[$n-1] + 1.0);
182             }
183             else
184             {
185 0         0 $y2[$n] = 0;
186             }
187              
188 1         3 for my $i (reverse 0 .. $n-1)
189             {
190 24         31 $y2[$i] = $y2[$i] * $y2[$i+1] + $u[$i];
191             }
192              
193 1         5 return @y2;
194             }
195              
196             =head3 Derivative1()
197              
198             A synonym for centraldiff().
199              
200             =cut
201              
202             #
203             # Alias Derivative1() to centraldiff(), and Derivative2() to
204             # seconddx(), preserving the old names. Not exporting the
205             # seconddx name now, as I'm not convinced it's a good name.
206             #
207             *Derivative1 = \¢raldiff;
208             *Derivative2 = \&seconddx;
209              
210             =head1 REFERENCES
211              
212             L
213              
214             L
215              
216             =head1 AUTHOR
217              
218             John A.R. Williams B
219              
220             John M. Gamble B (current maintainer)
221              
222             =cut
223              
224             1;