1 # Created by Octave 3.6.1, Sun Mar 11 22:05:14 2012 UTC <root@t61>
13 # name: <cell-element>
17 -- Function File: Y = Ci (Z)
18 Compute the cosine integral function defined by: Inf
24 See also: cosint, Si, sinint, expint, expint_Ei
30 # name: <cell-element>
34 Compute the cosine integral function defined by:
40 # name: <cell-element>
47 # name: <cell-element>
51 -- Function File: Y = Si (X)
52 Compute the sine integral defined by: x
61 # name: <cell-element>
65 Compute the sine integral defined by:
72 # name: <cell-element>
79 # name: <cell-element>
83 -- Function File: Y = cosint (Z)
84 Compute the cosine integral function defined by: Inf
86 cosint(x) = | cos(t)/t dt
90 See also: Ci, Si, sinint, expint, expint_Ei
96 # name: <cell-element>
100 Compute the cosine integral function defined by:
106 # name: <cell-element>
113 # name: <cell-element>
117 -- Function File: Y = dirac(X)
118 Compute the dirac delta function.
126 # name: <cell-element>
130 Compute the dirac delta function.
134 # name: <cell-element>
141 # name: <cell-element>
145 -- Function File: [K, E] = ellipke (M[,TOL])
146 Compute complete elliptic integral of first K(M) and second E(M).
148 M is either real array or scalar with 0 <= m <= 1
150 TOL will be ignored (MATLAB uses this to allow faster, less
151 accurate approximation)
153 Ref: Abramowitz, Milton and Stegun, Irene A. Handbook of
154 Mathematical Functions, Dover, 1965, Chapter 17.
162 # name: <cell-element>
166 Compute complete elliptic integral of first K(M) and second E(M).
170 # name: <cell-element>
177 # name: <cell-element>
181 -- Function File: erfcinv (X)
182 Compute the inverse complementary error function.
184 See also: erfc, erf, erfinv
190 # name: <cell-element>
194 Compute the inverse complementary error function.
198 # name: <cell-element>
205 # name: <cell-element>
209 -- Function File: Y = expint (X)
210 Compute the exponential integral, infinity
212 expint(x) = | exp(t)/t dt
216 See also: expint_E1, expint_Ei
222 # name: <cell-element>
226 Compute the exponential integral,
232 # name: <cell-element>
239 # name: <cell-element>
243 -- Function File: Y = expint_E1 (X)
244 Compute the exponential integral, infinity
246 expint(x) = | exp(t)/t dt
250 See also: expint, expint_Ei
256 # name: <cell-element>
260 Compute the exponential integral,
266 # name: <cell-element>
273 # name: <cell-element>
277 -- Function File: Y = expint_Ei (X)
278 Compute the exponential integral, infinity
280 expint_Ei(x) = - | exp(t)/t dt
284 See also: expint, expint_E1
290 # name: <cell-element>
294 Compute the exponential integral,
300 # name: <cell-element>
307 # name: <cell-element>
311 -- Function File: heaviside(X)
312 -- Function File: heaviside(X, ZERO_VALUE)
313 Compute the Heaviside step function.
315 The Heaviside function is defined as
317 Heaviside (X) = 1, X > 0
318 Heaviside (X) = 0, X < 0
320 The value of the Heaviside function at X = 0 is by default 0.5,
321 but can be changed via the optional second input argument.
329 # name: <cell-element>
333 Compute the Heaviside step function.
337 # name: <cell-element>
344 # name: <cell-element>
348 -- Function File: Y = laguerre (X,N)
349 -- Function File: [Y P]= laguerre (X,N)
350 Compute the value of the Laguerre polynomial of order N for each
357 # name: <cell-element>
361 Compute the value of the Laguerre polynomial of order N for each
367 # name: <cell-element>
374 # name: <cell-element>
378 -- Function File: X = lambertw (Z)
379 -- Function File: X = lambertw (Z, N)
380 Compute the Lambert W function of Z.
382 This function satisfies W(z).*exp(W(z)) = z, and can thus be used
383 to express solutions of transcendental equations involving
384 exponentials or logarithms.
386 N must be integer, and specifies the branch of W to be computed;
387 W(z) is a shorthand for W(0,z), the principal branch. Branches 0
388 and -1 are the only ones that can take on non-complex values.
390 If either N or Z are non-scalar, the function is mapped to each
391 element; both may be non-scalar provided their dimensions agree.
393 This implementation should return values within 2.5*eps of its
394 counterpart in Maple V, release 3 or later. Please report any
395 discrepancies to the author, Nici Schraudolph
396 <schraudo@inf.ethz.ch>.
398 For further details, see:
400 Corless, Gonnet, Hare, Jeffrey, and Knuth (1996), `On the Lambert
401 W Function', Advances in Computational Mathematics 5(4):329-359.
406 # name: <cell-element>
410 Compute the Lambert W function of Z.
414 # name: <cell-element>
421 # name: <cell-element>
425 LAPLACIAN Sparse Negative Laplacian in 1D, 2D, or 3D
427 [~,~,A]=LAPLACIAN(N) generates a sparse negative 3D Laplacian matrix
428 with Dirichlet boundary conditions, from a rectangular cuboid regular
429 grid with j x k x l interior grid points if N = [j k l], using the
430 standard 7-point finite-difference scheme, The grid size is always
431 one in all directions.
433 [~,~,A]=LAPLACIAN(N,B) specifies boundary conditions with a cell array
434 B. For example, B = {'DD' 'DN' 'P'} will Dirichlet boundary conditions
435 ('DD') in the x-direction, Dirichlet-Neumann conditions ('DN') in the
436 y-direction and period conditions ('P') in the z-direction. Possible
437 values for the elements of B are 'DD', 'DN', 'ND', 'NN' and 'P'.
439 LAMBDA = LAPLACIAN(N,B,M) or LAPLACIAN(N,M) outputs the m smallest
440 eigenvalues of the matrix, computed by an exact known formula, see
441 http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative
442 It will produce a warning if the mth eigenvalue is equal to the
443 (m+1)th eigenvalue. If m is absebt or zero, lambda will be empty.
445 [LAMBDA,V] = LAPLACIAN(N,B,M) also outputs orthonormal eigenvectors
446 associated with the corresponding m smallest eigenvalues.
448 [LAMBDA,V,A] = LAPLACIAN(N,B,M) produces a 2D or 1D negative
449 Laplacian matrix if the length of N and B are 2 or 1 respectively.
450 It uses the standard 5-point scheme for 2D, and 3-point scheme for 1D.
453 [lambda,V,A] = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
454 % Everything for 3D negative Laplacian with mixed boundary conditions.
455 laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
457 lambda = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
458 % computes the eigenvalues only
460 [~,V,~] = laplacian([200 200],{'DD' 'DN'},30);
461 % Eigenvectors of 2D negative Laplacian with mixed boundary conditions.
463 [~,~,A] = laplacian(200,{'DN'},30);
464 % 1D negative Laplacian matrix A with mixed boundary conditions.
466 % Example to test if outputs correct eigenvalues and vectors:
467 [lambda,V,A] = laplacian([13,10,6],{'DD' 'DN' 'P'},30);
468 [Veig D] = eig(full(A)); lambdaeig = diag(D(1:30,1:30));
469 max(abs(lambda-lambdaeig)) %checking eigenvalues
470 subspace(V,Veig(:,1:30)) %checking the invariant subspace
471 subspace(V(:,1),Veig(:,1)) %checking selected eigenvectors
472 subspace(V(:,29:30),Veig(:,29:30)) %a multiple eigenvalue
474 % Example showing equivalence between laplacian.m and built-in MATLAB
475 % DELSQ for the 2D case. The output of the last command shall be 0.
476 A1 = delsq(numgrid('S',32)); % input 'S' specifies square grid.
477 [~,~,A2] = laplacian([30,30]);
480 Class support for inputs:
481 N - row vector float double
483 M - scalar float double
485 Class support for outputs:
486 lambda and V - full float double, A - sparse float double.
488 Note: the actual numerical entries of A fit int8 format, but only
489 double data class is currently (2010) supported for sparse matrices.
491 This program is designed to efficiently compute eigenvalues,
492 eigenvectors, and the sparse matrix of the (1-3)D negative Laplacian
493 on a rectangular grid for Dirichlet, Neumann, and Periodic boundary
494 conditions using tensor sums of 1D Laplacians. For more information on
496 http://en.wikipedia.org/wiki/Kronecker_sum_of_discrete_Laplacians
497 For 2D case in MATLAB, see
498 http://www.mathworks.com/access/helpdesk/help/techdoc/ref/kron.html.
500 This code is also part of the BLOPEX package:
501 http://en.wikipedia.org/wiki/BLOPEX or directly
502 http://code.google.com/p/blopex/
506 # name: <cell-element>
510 LAPLACIAN Sparse Negative Laplacian in 1D, 2D, or 3D
515 # name: <cell-element>
522 # name: <cell-element>
526 -- Function File: [Y ALPHA] = multinom (X, N)
527 -- Function File: [Y ALPHA] = multinom (X, N,SORT)
528 Returns the terms (monomials) of the multinomial expansion of
531 (x1 + x2 + ... + xm)^N
533 X is a nT-by-m matrix where each column represents a different
534 variable, the output Y has the same format. The order of the
535 terms is inherited from multinom_exp and can be controlled through
536 the optional argument SORT and is passed to the function `sort'.
537 The exponents are returned in ALPHA.
539 See also: multinom_exp, multinom_coeff, sort
545 # name: <cell-element>
549 Returns the terms (monomials) of the multinomial expansion of degree n.
553 # name: <cell-element>
560 # name: <cell-element>
564 -- Function File: [C ALPHA] = multinom_coeff (M, N)
565 -- Function File: [C ALPHA] = multinom_coeff (M, N,ORDER)
566 Produces the coefficients of the multinomial expansion
568 (x1 + x2 + ... + xm).^n
570 For example, for m=3, n=3 the expansion is
573 = x1^3 + x2^3 + x3^3 +
574 + 3 x1^2 x2 + 3 x1^2 x3 + 3 x2^2 x1 + 3 x2^2 x3 +
575 + 3 x3^2 x1 + 3 x3^2 x2 + 6 x1 x2 x3
577 and the coefficients are [6 3 3 3 3 3 3 1 1 1].
579 The order of the coefficients is defined by the optinal argument
580 ORDER. It is passed to the function `multion_exp'. See the help
581 of that function for explanation. The multinomial coefficients
586 | | = ------------------------
587 | k | k(1)!k(2)! ... k(end)!
590 See also: multinom, multinom_exp
596 # name: <cell-element>
600 Produces the coefficients of the multinomial expansion
605 # name: <cell-element>
612 # name: <cell-element>
616 -- Function File: ALPHA = multinom_exp (M, N)
617 -- Function File: ALPHA = multinom_exp (M, N,SORT)
618 Returns the exponents of the terms in the multinomial expansion
620 (x1 + x2 + ... + xm).^N
622 For example, for m=2, n=3 the expansion has the terms
624 x1^3, x2^3, x1^2*x2, x1*x2^2
626 then `alpha = [3 0; 2 1; 1 2; 0 3]';
628 The optional argument SORT is passed to function `sort' to sort
629 the exponents by the maximum degree. The example above calling `
630 multinom(m,n,"ascend")' produces
632 `alpha = [2 1; 1 2; 3 0; 0 3]';
634 calling ` multinom(m,n,"descend")' produces
636 `alpha = [3 0; 0 3; 2 1; 1 2]';
638 See also: multinom, multinom_coeff, sort
644 # name: <cell-element>
648 Returns the exponents of the terms in the multinomial expansion
653 # name: <cell-element>
660 # name: <cell-element>
664 -- Function File: Y = psi (X)
665 Compute the psi function, for each value of X.
668 psi(x) = __ log(gamma(x))
671 See also: gamma, gammainc, gammaln
677 # name: <cell-element>
681 Compute the psi function, for each value of X.
685 # name: <cell-element>
692 # name: <cell-element>
696 -- Function File: Y = sinint (X)
697 Compute the sine integral function.
705 # name: <cell-element>
709 Compute the sine integral function.
713 # name: <cell-element>
720 # name: <cell-element>
724 -- Function File: Z = zeta (T)
725 Compute the Riemann's Zeta function.
733 # name: <cell-element>
737 Compute the Riemann's Zeta function.