Bessel functions of the second kind
Syntax
Definition
where 
A solution 
where 
and 
besselj.
Description
Y = bessely(nu,Z) computes Bessel functions of the second kind, 
Z. The order nu need not be an integer, but must be real. The argument Z can be complex. The result is real where Z is positive.
If nu and Z are arrays of the same size, the result is also that size. If either input is a scalar, it is expanded to the other input's size. If one input is a row vector and the other is a column vector, the result is a two-dimensional table of function values.
Y = bessely(nu,Z,1) computes bessely(nu,Z).*exp(-abs(imag(Z))).
[Y,ierr] = bessely(nu,Z) also returns completion flags in an array the same size as Y.
Remarks
The Bessel functions are related to the Hankel functions, also called Bessel functions of the third kind,
where 
besselh, 
besselj, and 
bessely. The Hankel functions also form a fundamental set of solutions to Bessel's equation (see besselh).
Examples
format long
z = (0:0.2:1)';
bessely(1,z)
ans =
-Inf
-3.32382498811185
-1.78087204427005
-1.26039134717739
-0.97814417668336
-0.78121282130029
Example 2. bessely(3:9,(0:.2:10)') generates the entire table on page 399 of [1] Abramowitz and Stegun, Handbook of Mathematical Functions.
Algorithm
The bessely function uses a Fortran MEX-file to call a library developed by D. E Amos [3] [4].
Syntax
Definition
is the gamma function.Description
B = beta(Z,W) computes the beta function for corresponding elements of arrays Z and W. The arrays must be real and nonnegative. They must be the same size, or either can be scalar.
Examples
In this example, which uses integer arguments,
is the ratio of fairly small integers, and the rational format is able to recover the exact result.
Algorithm
Syntax
Definition
The incomplete beta function is
where 
is the gamma function.Description
I = betainc(X,Z,W) computes the incomplete beta function for corresponding elements of the arrays X, Z and W. The elements of X must be in the closed interval 
Z and W must be nonnegative and real. All arrays must be the same size, or any of them can be scalar.
Examples
format long
betainc(.5,(0:10)',3)
ans =
1.00000000000000
0.87500000000000
0.68750000000000
0.50000000000000
0.34375000000000
0.22656250000000
0.14453125000000
0.08984375000000
0.05468750000000
0.03271484375000
0.01928710937500
Syntax
Description
L = betaln(Z,W) computes the natural logarithm of the beta function log(beta(Z,W)), for corresponding elements of arrays Z and W, without computing beta(Z,W). Since the beta function can range over very large or very small values, its logarithm is sometimes more useful.
Z and W must be real and nonnegative. They must be the same size, or either can be scalar.
Examples
-708.8616 is slightly less than log(realmin). Computing beta(x,x) directly would underflow (or be denormal).
Algorithm
Syntax
x = bicg(A,b)
bicg(A,b,tol)
bicg(A,b,tol,maxit)
bicg(A,b,tol,maxit,M)
bicg(A,b,tol,maxit,M1,M2)
bicg(A,b,tol,maxit,M1,M2,x0)
bicg(afun,b,tol,maxit,mfun1,mfun2,x0,p1,p2,...)
[x,flag] = bicg(A,b,...)
[x,flag,relres] = bicg(A,b,...)
[x,flag,relres,iter] = bicg(A,b,...)
[x,flag,relres,iter,resvec] = bicg(A,b,...)
Description
x = bicg(A,b) attempts to solve the system of linear equations A*x = b for x. The n-by-n coefficient matrix A must be square and should be large and sparse. The column vector b must have length n. A can be a function afun such that afun(x) returns A*x and afun(x,'transp') returns A'*x.
If bicg converges, it displays a message to that effect. If bicg fails to converge after the maximum number of iterations or halts for any reason, it prints a warning message that includes the relative residual norm(b-A*x)/norm(b) and the iteration number at which the method stopped or failed.
bicg(A,b,tol) specifies the tolerance of the method. If tol is [], then bicg uses the default, 1e-6.
bicg(A,b,tol,maxit) specifies the maximum number of iterations. If maxit is [], then bicg uses the default, min(n,20).
bicg(A,b,tol,maxit,M) and bicg(A,b,tol,maxit,M1,M2) use the preconditioner M or M = M1*M2 and effectively solve the system inv(M)*A*x = inv(M)*b for x. If M is [] then bicg applies no preconditioner. M can be a function mfun such that mfun(x) returns M\x and mfun(x,'transp') returns M'\x.
bicg(A,b,tol,maxit,M1,M2,x0) specifies the initial guess. If x0 is [], then bicg uses the default, an all-zero vector.
bicg(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,...) passes parameters p1,p2,... to functions afun(x,p1,p2,...) and afun(x,p1,p2,...,'transp'), and similarly to the preconditioner functions m1fun and m2fun.
[x,flag] = bicg(A,b,...) also returns a convergence flag.
Whenever flag is not 0, the solution x returned is that with minimal norm residual computed over all the iterations. No messages are displayed if the flag output is specified.
[x,flag,relres] = bicg(A,b,...) also returns the relative residual norm(b-A*x)/norm(b). If flag is 0, relres <= tol.
[x,flag,relres,iter] = bicg(A,b,...) also returns the iteration number at which x was computed, where 0 <= iter <= maxit.
[x,flag,relres,iter,resvec] = bicg(A,b,...) also returns a vector of the residual norms at each iteration including norm(b-A*x0).
Examples
n = 100;
on = ones(n,1);
A = spdiags([-2*on 4*on -on],-1:1,n,n);
b = sum(A,2);
tol = 1e-8;
maxit = 15;
M1 = spdiags([on/(-2) on],-1:0,n,n);
M2 = spdiags([4*on -on],0:1,n,n);
x = bicg(A,b,tol,maxit,M1,M2,[]);
Alternatively, use this matrix-vector product function
function y = afun(x,n,transp_flag)
if (nargin > 2) & strcmp(transp_flag,'transp')
y = 4 * x;
y(1:n-1) = y(1:n-1) - 2 * x(2:n);
y(2:n) = y(2:n) - x(1:n-1);
else
y = 4 * x;
y(2:n) = y(2:n) - 2 * x(1:n-1);
y(1:n-1) = y(1:n-1) - x(2:n);
end
Example 2. This examples demonstrates the use of a preconditioner. Start with A = west0479, a real 479-by-479 sparse matrix, and define b so that the true solution is a vector of all ones.
You can accurately solve A*x = b using backslash since A is not so large.
Now try to solve A*x = b with bicg.
The value of flag indicates that bicg iterated the default 20 times without converging. The value of iter shows that the method behaved so badly that the initial all-zero guess was better than all the subsequent iterates. The value of relres supports this: relres = norm(b-A*x)/norm(b) = norm(b)/norm(b) = 1. You can confirm that the unpreconditioned method oscillates rather wildly by plotting the relative residuals at each iteration.
Now, try an incomplete LU factorization with a drop tolerance of 1e-5 for the preconditioner.
[L1,U1] = luinc(A,1e-5);
Warning: Incomplete upper triangular factor has 1 zero diagonal.
It cannot be used as a preconditioner for an iterative
method.
nnz(A), nnz(L1), nnz(U1)
ans =
1887
ans =
5562
ans =
4320
The zero on the main diagonal of the upper triangular U1 indicates that U1 is singular. If you try to use it as a preconditioner,
the method fails in the very first iteration when it tries to solve a system of equations involving the singular U1 using backslash. bicg is forced to return the initial estimate since no other iterates were produced.
Try again with a slightly less sparse preconditioner.
This time U2 is nonsingular and may be an appropriate preconditioner.
and bicg converges to within the desired tolerance at iteration number 8. Decreasing the value of the drop tolerance increases the fill-in of the incomplete factors but also increases the accuracy of the approximation to the original matrix. Thus, the preconditioned system becomes closer to inv(U)*inv(L)*L*U*x = inv(U)*inv(L)*b, where L and U are the true LU factors, and closer to being solved within a single iteration.
The next graph shows the progress of bicg using six different incomplete LU factors as preconditioners. Each line in the graph is labeled with the drop tolerance of the preconditioner used in bicg.
No comments:
Post a Comment