|
5827
|
1 ## Copyright (C) 2000 Kai Habel |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
8 ## any later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301, USA. |
|
|
19 |
|
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {@var{L} =} legendre (@var{n}, @var{X}) |
|
|
22 ## |
|
|
23 ## Legendre Function of degree n and order m |
|
|
24 ## where all values for m = 0..@var{n} are returned. |
|
|
25 ## @var{n} must be a scalar in the range [0..255]. |
|
|
26 ## The return value has one dimension more than @var{x}. |
|
|
27 ## |
|
|
28 ## @example |
|
|
29 ## The Legendre Function of degree n and order m |
|
|
30 ## |
|
|
31 ## @group |
|
|
32 ## m m 2 m/2 d^m |
|
|
33 ## P(x) = (-1) * (1-x ) * ---- P (x) |
|
|
34 ## n dx^m n |
|
|
35 ## @end group |
|
|
36 ## |
|
|
37 ## with: |
|
|
38 ## Legendre Polynom of degree n |
|
|
39 ## |
|
|
40 ## @group |
|
|
41 ## 1 d^n 2 n |
|
|
42 ## P (x) = ------ [----(x - 1) ] |
|
|
43 ## n 2^n n! dx^n |
|
|
44 ## @end group |
|
|
45 ## |
|
|
46 ## legendre(3,[-1.0 -0.9 -0.8]) returns the matrix |
|
|
47 ## |
|
|
48 ## @group |
|
|
49 ## x | -1.0 | -0.9 | -0.8 |
|
|
50 ## ------------------------------------ |
|
|
51 ## m=0 | -1.00000 | -0.47250 | -0.08000 |
|
|
52 ## m=1 | 0.00000 | -1.99420 | -1.98000 |
|
|
53 ## m=2 | 0.00000 | -2.56500 | -4.32000 |
|
|
54 ## m=3 | 0.00000 | -1.24229 | -3.24000 |
|
|
55 ## @end group |
|
|
56 ## @end example |
|
|
57 ## @end deftypefn |
|
|
58 |
|
|
59 ## FIXME Add Schmidt semi-normalized and fully normalized legendre functions |
|
|
60 |
|
|
61 ## Author: Kai Habel <kai.habel@gmx.de> |
|
|
62 |
|
|
63 function L = legendre (n, x) |
|
|
64 |
|
|
65 warning ("legendre is unstable for higher orders"); |
|
|
66 |
|
|
67 if (nargin != 2) |
|
|
68 print_usage (); |
|
|
69 endif |
|
|
70 |
|
|
71 if (! isscalar (n) || n < 0 || n > 255 || n != fix (n)) |
|
|
72 error ("n must be a integer between 0 and 255]"); |
|
|
73 endif |
|
|
74 |
|
|
75 if (! isvector (x) || any (x < -1 || x > 1)) |
|
|
76 error ("x must be vector in range -1 <= x <= 1"); |
|
|
77 endif |
|
|
78 |
|
|
79 if (n == 0) |
|
|
80 L = ones (size (x)); |
|
|
81 elseif (n == 1) |
|
|
82 L = [x; -sqrt(1 - x .^ 2)]; |
|
|
83 else |
|
|
84 i = 0:n; |
|
|
85 a = (-1) .^ i .* bincoeff (n, i); |
|
|
86 p = [a; zeros(size (a))]; |
|
|
87 p = p(:); |
|
|
88 p(length (p)) = []; |
|
|
89 #p contains the polynom (x^2-1)^n |
|
|
90 |
|
|
91 #now create a vector with 1/(2.^n*n!)*(d/dx).^n |
|
|
92 d = [((n + rem(n, 2)):-1:(rem (n, 2) + 1)); 2 * ones(fix (n / 2), n)]; |
|
|
93 d = cumsum (d); |
|
|
94 d = fliplr (prod (d')); |
|
|
95 d = [d; zeros(1, length (d))]; |
|
|
96 d = d(1:n + 1) ./ (2 ^ n *prod (1:n)); |
|
|
97 |
|
|
98 Lp = d' .* p(1:length (d)); |
|
|
99 #Lp contains the Legendre Polynom of degree n |
|
|
100 |
|
|
101 # now create a polynom matrix with d/dx^m for m=0..n |
|
|
102 d2 = flipud (triu (ones (n))); |
|
|
103 d2 = cumsum (d2); |
|
|
104 d2 = fliplr (cumprod (flipud (d2))); |
|
|
105 d3 = fliplr (triu (ones (n + 1))); |
|
|
106 d3(2:n + 1, 1:n) = d2; |
|
|
107 |
|
|
108 # multiply for each m (d/dx)^m with Lp(n,x) |
|
|
109 # and evaluate at x |
|
|
110 Y = zeros(n + 1, length (x)); |
|
|
111 [dr, dc] = size (d3); |
|
|
112 for m = 0:dr - 1 |
|
|
113 Y(m + 1, :) = polyval (d3(m + 1, 1:(dc - m)) .* Lp(1:(dc - m))', x)(:)'; |
|
|
114 endfor |
|
|
115 |
|
|
116 # calculate (-1)^m*(1-x^2)^(m/2) for m=0..n at x |
|
|
117 # and multiply with (d/dx)^m(Pnx) |
|
|
118 l = length (x); |
|
|
119 X = kron ((1 - x(:) .^ 2)', ones (n + 1, 1)); |
|
|
120 M = kron ((0:n)', ones (1, l)); |
|
|
121 L = X .^ (M / 2) .* (-1) .^ M .* Y; |
|
|
122 endif |
|
|
123 endfunction |
|
|
124 |
|
|
125 %!test |
|
|
126 %! result=legendre(3,[-1.0 -0.9 -0.8]); |
|
|
127 %! expected = [ |
|
|
128 %! -1.00000 -0.47250 -0.08000 |
|
|
129 %! 0.00000 -1.99420 -1.98000 |
|
|
130 %! 0.00000 -2.56500 -4.32000 |
|
|
131 %! 0.00000 -1.24229 -3.24000 |
|
|
132 %! ]; |
|
|
133 %! assert(result,expected,1e-5); |
