Mercurial > matrix-functions
view toolbox/wathen.m @ 2:c124219d7bfa draft
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Thu, 07 May 2015 18:36:24 +0200 |
parents | 8f23314345f4 |
children |
line wrap: on
line source
function A = wathen(nx, ny, k) %WATHEN Wathen matrix - a finite element matrix (sparse, random entries). % A = WATHEN(NX,NY) is a sparse random N-by-N finite element matrix % where N = 3*NX*NY + 2*NX + 2*NY + 1. % A is precisely the `consistent mass matrix' for a regular NX-by-NY % grid of 8-node (serendipity) elements in 2 space dimensions. % A is symmetric positive definite for any (positive) values of % the `density', RHO(NX,NY), which is chosen randomly in this routine. % In particular, if D = DIAG(DIAG(A)), then % 0.25 <= EIG(INV(D)*A) <= 4.5 % for any positive integers NX and NY and any densities RHO(NX,NY). % This diagonally scaled matrix is returned by WATHEN(NX,NY,1). % Reference: % A.J. Wathen, Realistic eigenvalue bounds for the Galerkin % mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449-457. if nargin < 2, error('Two dimensioning arguments must be specified.'), end if nargin < 3, k = 0; end e1 = [6 -6 2 -8;-6 32 -6 20;2 -6 6 -6;-8 20 -6 32]; e2 = [3 -8 2 -6;-8 16 -8 20;2 -8 3 -8;-6 20 -8 16]; e = [e1 e2; e2' e1]/45; n = 3*nx*ny+2*nx+2*ny+1; A = sparse(n,n); RHO = 100*rand(nx,ny); for j=1:ny for i=1:nx nn(1) = 3*j*nx+2*i+2*j+1; nn(2) = nn(1)-1; nn(3) = nn(2)-1; nn(4) = (3*j-1)*nx+2*j+i-1; nn(5) = 3*(j-1)*nx+2*i+2*j-3; nn(6) = nn(5)+1; nn(7) = nn(6)+1; nn(8) = nn(4)+1; em = e*RHO(i,j); for krow=1:8 for kcol=1:8 A(nn(krow),nn(kcol)) = A(nn(krow),nn(kcol))+em(krow,kcol); end end end end if k == 1 A = diag(diag(A)) \ A; end