% % Problem 1 % n=3; B=4*diag(ones(1,n))-diag(ones(1,n-1),1)-diag(ones(1,n-1),-1) O=zeros(n,n) mI=-eye(n,n) % Actually, it is not easy to build a matrix out of block matrices % Or does anybody know how. (Compare Maple's band-command) % Anyway, it is more efficient to form it diagonal by diagonal. % Let m be the nr. of subintervals % n=m-1 % nr. of interior grid point in each axis % N=N^2 % Nr. of inner nodes on the 2-dim. grid = nr. of unknowns m=4;n=m-1;N=n^2 maind= % help ones neardiag= ... fardiag= ... diag(maind) diag(neardiag,1) diag(neardiag,-1) diag(fardiag,n) diag(fardiag,-n) diag(maind)+... diag(neardiag,1) +... diag(neardiag,-1) + ... diag(fardiag,n) +... diag(fardiag,-n) % It's good, except you need to have 0 at each nth position of neardiag. % help rem % 1:N-1 % gives [1,2,...,N-1] % use == or ~= % Note: scalar extension and comparison componentwise neardiag=... A=... diag(maind)+... diag(neardiag,1) +... diag(neardiag,-1) + ... diag(fardiag,n) +... diag(fardiag,-n) %%%%%%%%%%%%%%%%%% % Sparse matrices: %help sparse %help spdiags % The easiest way to convert the matrix in sparse storage is %A=sparse(A) % This is not an economic way to build a sparse matrix though. % (Often the full matrix doesn't fit into memory) % % spdiags offers a way to build banded sparse matices storing just the % diagonals. % m=4;n=m-1;N=n^2 maind=4*ones(N,1) fardiag=-ones(N,1) neardiag=rem(1:N,n)~=0 [fardiag,neardiag',maind,neardiag',fardiag] % This looks good but isn't quite. There's a little trick about spdiags % not explained in help spdiags. It uses the data % starting at the top on the LHS (below) the maindiag but % starting at the bottom on the RHS (above) the maindiag. % (This doesn't show on constant diagonals.) % Thus we need to append 1 at the head and drop one extra element at the % tail of the upper neardiag. uneardiag=[1,neardiag];uneardiag=uneardiag(1:N); [fardiag,neardiag',maind,uneardiag',fardiag] A=spdiags([fardiag,neardiag',maind,uneardiag',fardiag],[-n,-1,0,1,n],N,N) full(A) % If you want to see a part of it type % AF=full(A); % AF(1:5,1:5) % for instance %Now, let's write it as a function: %Save it in the file lapm.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Problem nr. 2 %%%%%%%%%%%%%%% % left=[20,40,60,80];top=[50,20,10,0];bot=zeros(1,4),right=left m=5;n=m-1; %First build B as a 4x4 matrix. B(:) makes it a vector (columnwise). %Thus you need to transpose B before "raveling" it into a long vector. A=lapm(m); % if you did'n write it as a function, cut and paste the commands from above. Tvek= ... Use \ (backslash) T=reshape(Tvek,4,4)' % reshape works columnwise, thus you need to transpose T=flipud(T) % Flip to fit the coordinate-axes T=[fliplr(left)' T fliplr(right)'] % Vertical boundaries need flips as well NW=0.5*(left(n)+top(1));NE=0.5*(top(n)+right(n)); SE=0.5*(bot(n)+right(1));SW=0.5*(bot(1)+left(1)); % When you are done, check it! for i=2:n; for j=2:n ; tt(i,j)= 4*T(i,j)-T(i-1,j)-T(i+1,j)-T(i,j-1)-T(i,j+1); end; end tt % There's a function for this: del2(T)