www.gusucode.com > elmat工具箱matlab源码程序 > elmat/private/wathen.m
function A = wathen(nx, ny, k, classname) %WATHEN Wathen matrix (sparse). % A = GALLERY('WATHEN',NX,NY) is a sparse matrix with random % elements. It is an 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 elements in two space dimensions. A is symmetric % positive definite for any (positive) values of the "density", % RHO(NX,NY), which is chosen randomly. 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). % % B = GALLERY('WATHEN',NX,NY,1) returns B = DIAG(DIAG(A)) \ A, % which is A row-scaled to have unit diagonal. % Reference: % A. J. Wathen, Realistic eigenvalue bounds for the Galerkin % mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449-457. % % Nicholas J. Higham % Revised by Tim Davis. % Copyright 1984-2007 The MathWorks, Inc. if isempty(ny) error(message('MATLAB:wathen:NotEnoughInputs')) end if isempty(k), 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; ntriplets = nx*ny*64 ; I = zeros (ntriplets, 1) ; J = zeros (ntriplets, 1) ; X = zeros (ntriplets, 1) ; ntriplets = 0 ; RHO = 100*rand(nx,ny); nn = zeros(8,1); 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 ntriplets = ntriplets + 1 ; I (ntriplets) = nn (krow) ; J (ntriplets) = nn (kcol) ; X (ntriplets) = em (krow,kcol) ; end end end end A = sparse (I,J,X,n,n) ; if k == 1 A = diag(diag(A)) \ A; end