A First Course in Numerical Methods

1. function [A,p] = house(A)
   
2. % function [A,p] = house(A)
   
3. %
   
4. % perform QR decomposition using Householder reflections
   
5. % Transformations are of the form P_k = I - 2u_k(u_k^T), so
   
6. % store effecting vector u_k in p(k) + A(k+1:m,k). Assume m > n.
   
7. 
   
8. [m,n]=size(A); p = zeros(1,n);
   
9. for k = 1:n
   
10.   % define u of length = m-k+1
    
11.   z = A(k:m,k);
    
12.   e1 = [1; zeros(m-k,1)];
    
13.   u = z+sign(z(1))*norm(z)*e1; u = u/norm(u);
    
14.   % update nonzero part of A by I-2uu^T
    
15.   A(k:m,k:n) = A(k:m,k:n)-2*u*(u'*A(k:m,k:n));
    
16.   % store u
    
17.   p(k) = u(1);
    
18.   A(k+1:m,k) = u(2:m-k+1);
    
19. end

复制代码

1. function coefs = lsfit (t, b, n)
   
2. %
   
3. % function coefs = lsfit (t, b, n)
   
4. %
   
5. % Construct coefficients of the polynomial of
   
6. % degree at most n-1 that best fits data (t,b)
   
7. 
   
8. t = t(:); b = b(:); % make sure t and b are column vectors
   
9. m = length(t);
   
10. 
    
11. % long and skinny A
    
12. A = ones(m,n);
    
13. for j=1:n-1
    
14.   A(:,j+1) = A(:,j).*t;
    
15. end
    
16. 
    
17. % normal equations and solution
    
18. B = A'*A; y = A'*b;
    
19. coefs = B \ y;

复制代码

1. function [x,nr] = lsol(b,AA,p)
   
2. %
   
3. % function [x,nr] = lsol(b,AA,p)
   
4. %
   
5. % Given the output AA and p of house, which contain Q and R of matrix A,
   
6. % and given a right-hand-side vector b,  solve min || b - Ax ||.
   
7. % Return also the norm of the residual, nr = ||r|| = ||b - Ax||.
   
8. 
   
9. y = b(:); [m,n] = size(AA);
   
10. % transform b
    
11. for k=1:n
    
12.   u = [p(k);AA(k+1:m,k)];
    
13.   y(k:m) = y(k:m) - 2*u *(u'*y(k:m));
    
14. end
    
15. % form upper triangular R and solve
    
16. R = triu(AA(1:n,:));
    
17. x = R \ y(1:n); nr = norm(y(n+1:m));

复制代码

1. function A=kron_conv_diff(beta,gamma,N)
   
2. %
   
3. % function A=kron_conv_diff(beta,gamma,N)
   
4. 
   
5. ee=ones(N,1);
   
6. a=4; b=-1-gamma; c=-1-beta; d=-1+beta; e=-1+gamma;
   
7. t1=spdiags([c*ee,a*ee,d*ee],-1:1,N,N);
   
8. t2=spdiags([b*ee,zeros(N,1),e*ee],-1:1,N,N);
   
9. A=kron(speye(N),t1)+kron(t2,speye(N));

复制代码

1. function [Q,T] = lanczos(A,Q,k)
   
2. %
   
3. % function [Q,T] = lanczos(A,Q,k)
   
4. %
   
5. 
   
6. % preallocate for speed
   
7. alpha=zeros(k,1); beta=zeros(k,1);
   
8. 
   
9. Q(:,1) = Q(:,1)/norm(Q(:,1));
   
10. beta(1,1)=0;
    
11. 
    
12. for j=1:k
    
13.   w=A*Q(:,j);
    
14.   if j>1
    
15.     w=A*Q(:,j)-beta(j,1)*Q(:,j-1);
    
16.   end
    
17.   alpha(j,1)=Q(:,j)'*w;
    
18.   w=w-alpha(j,1)*Q(:,j);
    
19.   beta(j+1,1)=norm(w);
    
20. 
    
21.   if abs(beta(j+1,1))<1e-10
    
22.     disp('Zero beta --- returning.');
    
23.     T=spdiags([beta(2:j+1) alpha(1:j) beta(1:j)],-1:1,j+1,j);
    
24.     return
    
25.   end
    
26.   Q(:,j+1)=w/beta(j+1,1);
    
27. end
    
28. T=spdiags([beta(2:end) alpha beta(1:end-1)],-1:1,k+1,k);

复制代码

1. function [x,iter,res] = pcgmg (A,b,x0,tol)
   
2. %
   
3. % function [x,iter,res] = pcgmg (A,b,x0,tol)
   
4. %
   
5. % CG with V-cycle preconditioning; assume x0=0
   
6. 
   
7. N = sqrt(length(b));
   
8. flevel = log2(N+1);
   
9. tol2 = tol^2;
   
10. x = x0;
    
11. r = b - A*x0;
    
12. h = poismg(A,r,x0,flevel);
    
13. d = r'*h; bb = d;
    
14. p = h;
    
15. 
    
16. iter = 0;
    
17. while d > tol2 * bb
    
18.   do = d;
    
19.   s = A*p;
    
20.   alfa = d / (p'*s);
    
21.   x = x + alfa*p;
    
22.   r = r - alfa*s;
    
23.   h = poismg(A,r,x0,flevel);
    
24.   d = r'*h;
    
25.   p = h + d/do*p;
    
26.   iter = iter + 1;
    
27.   res(iter) = norm(r);
    
28. end

复制代码

1. function A=kron_conv_diff(beta,gamma,N)
   
2. %
   
3. % function A=kron_conv_diff(beta,gamma,N)
   
4. 
   
5. ee=ones(N,1);
   
6. a=4; b=-1-gamma; c=-1-beta; d=-1+beta; e=-1+gamma;
   
7. t1=spdiags([c*ee,a*ee,d*ee],-1:1,N,N);
   
8. t2=spdiags([b*ee,zeros(N,1),e*ee],-1:1,N,N);
   
9. A=kron(speye(N),t1)+kron(t2,speye(N));

复制代码

function [A,p] = house(A)
% function [A,p] = house(A)
%
% perform QR decomposition using Householder reflections
% Transformations are of the form P_k = I - 2u_k(u_k^T), so
% store effecting vector u_k in p(k) + A(k+1:m,k). Assume m > n.

[m,n]=size(A); p = zeros(1,n);
for k = 1:n
  % define u of length = m-k+1
  z = A(k:m,k);
  e1 = [1; zeros(m-k,1)];
  u = z+sign(z(1))*norm(z)*e1; u = u/norm(u);
  % update nonzero part of A by I-2uu^T
  A(k:m,k:n) = A(k:m,k:n)-2*u*(u'*A(k:m,k:n));
  % store u
  p(k) = u(1);
  A(k+1:m,k) = u(2:m-k+1);
end

谢谢分享！~~！

谢谢分享～～！