Scientific Computing with Case Studies

Written for advanced undergraduate and early graduate courses in numerical analysis and scientific computing, this book provides a practical guide to the numerical solution of linear and nonlinear equations, differential equations, optimization problems, and eigenvalue problems. The book treats standard problems and introduces important variants such as sparse systems, differential-algebraic equations, constrained optimization, Monte Carlo simulations, and parametric studies.

In addition, a supplemental set of MATLAB code files is available for download. Optimization Toolbox is also discussed.

Retrieve Companion Software

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %                                                        %
   
3. %  A solution to CSE For Your Homework Project 1         %
   
4. %  Image Deblurring: I Can See Clearly Now               %
   
5. %  Dianne P. O'Leary and James G. Nagy 09/02             %
   
6. %                                                        %
   
7. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
8. 
   
9. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
10. %                                                        %
    
11. %  Read in a blurred image G and matrices A and B.       %
    
12. %  Try to reconstruct F, where g is approximately        %
    
13. %  kron(A,B) * reshape(F,m*n,1) and g = reshape(G,m*n,1) %
    
14. %                                                        %
    
15. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
16. 
    
17. load proj1data.mat
    
18. 
    
19. imagesc(G)
    
20. colormap(gray)
    
21. figure(2)
    
22. 
    
23. [m, n] = size(G);
    
24. 
    
25. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
26. %                                                        %
    
27. %  Compute the singular value decomposition of A and B.  %
    
28. %                                                        %
    
29. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
30. 
    
31. [Ua, Sa, Va] = svd(A);
    
32. [Ub, Sb, Vb] = svd(B);
    
33. 
    
34. Ghat = Ub'*G*Ua;
    
35. 
    
36. S = diag(Sb)*(diag(Sa))';
    
37. 
    
38. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
39. %                                                        %
    
40. %  Use Tikhonov regularization, with parameters          %
    
41. %  specified in lamtest, to reconstruct the image.       %
    
42. %                                                        %
    
43. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
44. 
    
45. lamtest = [.05 .01 .005 .0025 .0015 .001 .00095 .00016681 .00005];
    
46. 
    
47. for lam=lamtest
    
48.   Fhat = (S.*Ghat) ./ (S.*S+lam^2);
    
49.   F = Vb*Fhat*Va';
    
50.   imagesc(F)
    
51.   colormap(gray)
    
52.   disp(sprintf('Tikhonov lambda= %f',lam))
    
53.   disp('Strike any key to continue.')
    
54.   drawnow
    
55.   pause
    
56. end
    
57. 
    
58. 
    
59. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
60. %                                                        %
    
61. %  Use TSVD regularization to reconstruct the image.     %
    
62. %                                                        %
    
63. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
64. 
    
65. % Store the entries of S in a linear array s, then sort this
    
66. % set in descending order and construct the permutation
    
67. % that locates the ordered set in the unsorted list.
    
68. 
    
69. s = reshape(S,n^2,1);
    
70. [ss,prm] = sort(s);
    
71. prm = flipud(prm);
    
72. ss = flipud(ss);
    
73. ls = length(s);
    
74. iprm = zeros(ls,1);
    
75. iprm(prm) = [1:ls]';
    
76. 
    
77. 
    
78. % Using the sorted set and permutation constructed above,
    
79. % remove all but the p largest singular values and then
    
80. % compute the truncated expanded solution.  Test using a
    
81. % set of values of p.
    
82. 
    
83. pset = [100, 250, 500, 1000, 1500, 2000, 2500, 4000, 6000, 8000];
    
84. 
    
85. for i=1:length(pset),
    
86.   p = pset(i);
    
87.   ssnew = [ss(1:p); zeros(length(ss)-p,1)];
    
88.   Snew = reshape(ssnew(iprm),n,n);
    
89.   Fhat = Ghat ./ S;
    
90.   Fnew = Fhat .* (Snew>0);
    
91.   F = Vb*Fnew*Va';
    
92.   imagesc(F)
    
93.   disp(sprintf('TSVD p= %d ', p))
    
94.   disp('Strike any key to continue.')
    
95.   drawnow
    
96.   pause;
    
97. end
    
98. 
    
99. 
    
100. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
101. %                                                        %
     
102. %  Here's an alternative approach to TSVD regularization,%
     
103. %  which determines p implicitly by throwing out all     %
     
104. %  singular values less than some value, here            %
     
105. %  determined by the parameters in lamtest.              %
     
106. %                                                        %
     
107. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
108. 
     
109. Fhat = Ghat ./ S;
     
110. for lam=lamtest
     
111.   ind = S > lam;
     
112.   Fnew = Fhat .* ind;
     
113.   F = Vb*Fnew*Va';
     
114.   imagesc(F)
     
115.   disp(sprintf('TSVD lambda = %f, p = %d',lam,sum(sum(ind))))
     
116.   disp('Strike any key to continue.')
     
117.   drawnow
     
118.   pause
     
119. end

复制代码

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %                                                      %
   
3. % This sample code demonstrates how to manipulate      %
   
4. % Kronecker products of matrices without using a lot   %
   
5. % of storage.                                          %
   
6. %                                                      %
   
7. % It solves the linear system  K f = g                 %
   
8. % where K = kron(A,B)                                  %
   
9. % using the singular value decompositions of A and B.  %
   
10. %                                                      %
    
11. % Notation:  F is a square array.  f contains these    %
    
12. %                same entries, stacked by column.      %
    
13. %            G is a square array.  g contains these    %
    
14. %                same entries, stacked by column.      %
    
15. %                                                      %
    
16. % projdemo.m     Dianne P O'Leary 09/02                %
    
17. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
18. 
    
19. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
20. %                                                      %
    
21. % Initialize the data.  The problem is small enough    %
    
22. % so that you can print everything out and look at it. %
    
23. %                                                      %
    
24. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
25. 
    
26. n = 3;
    
27. F = rand(n,n);
    
28. A = rand(n,n);
    
29. B = rand(n,n) + eye(n);
    
30. 
    
31. 
    
32. f = reshape(F,n*n,1);
    
33. K = kron(A,B);
    
34. g = K * f;
    
35. 
    
36. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
37. %                                                      %
    
38. % We could solve the problem by saying fcomp = K \ g,  %
    
39. % but instead we will take advantage of the Kronecker  %
    
40. % structure to avoid ever needing K or its factors.    %
    
41. %                                                      %
    
42. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
43. 
    
44. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
45. %                                                      %
    
46. % Here is another way to compute the right-hand side.  %
    
47. %                                                      %
    
48. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
49. 
    
50. G = B * F * A';
    
51. 
    
52. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
53. %                                                      %
    
54. % Notice that it gives the same data, so the           %
    
55. % following norm prints as a very small number.        %
    
56. %                                                      %
    
57. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
58. 
    
59. Gg = reshape(g,n,n);
    
60. disp('norm(G - Gg)')
    
61. disp(norm(G - Gg, 'fro'))
    
62. 
    
63. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
64. %                                                      %
    
65. % Let's solve the problem using SVD's of A and B,      %
    
66. % as a model for understanding the data manipulation   %
    
67. % in your homework.                                    %
    
68. %                                                      %
    
69. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
70. 
    
71. 
    
72. [Ua, Sa, Va] = svd(A);
    
73. [Ub, Sb, Vb] = svd(B);
    
74. 
    
75. Ghat = Ub'*G*Ua;
    
76. 
    
77. S = diag(Sb)*(diag(Sa))';
    
78. 
    
79. Fhat = Ghat ./ S;
    
80. Fc = Vb*Fhat*Va';
    
81. 
    
82. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
83. %                                                      %
    
84. % Notice that it gives almost the exact answer, so the %
    
85. % following norm prints as a very small number.        %
    
86. %                                                      %
    
87. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
88. 
    
89. disp('norm(F - Fc)')
    
90. disp(norm(F - Fc, 'fro'))
    
91. 
    
92. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
93. %                                                      %
    
94. % For completeness, we verify that                     %
    
95. % kron(Ua,Ub) S kron(Va,Vb) is the SVD of K.           %
    
96. %                                                      %
    
97. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
98. 
    
99. 
    
100. disp('norm(K-svd(K))')
     
101. disp(norm(K - kron(Ua,Ub)*diag(reshape(S,n*n,1))*kron(Va,Vb)', 'fro'))
     
102. disp('norm(UU^{T}-I)')
     
103. disp(norm(kron(Ua,Ub)*kron(Ua,Ub)' - eye(n*n),'fro'))
     
104. disp('norm(VV^{T}-I)')
     
105. disp(norm(kron(Va,Vb)*kron(Va,Vb)' - eye(n*n),'fro'))

复制代码

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %  Solution to Problem 3 in the homework assignment,
   
3. %  Updating and Downdating Matrix Factorizations:
   
4. %     A Change in Plans
   
5. %
   
6. %  We solve two truss problems, the first through matrix
   
7. %  factorization and the second using the Sherman-Morrison-Woodbury
   
8. %  formula.
   
9. %
   
10. %  First truss:  All members are equal length.  
    
11. %                Assume a load of 50 at C.
    
12. %
    
13. %                  B           D
    
14. %                 / \         / \
    
15. %                /   \       /   \
    
16. %               /     \     /     \
    
17. %              /       \   /       \
    
18. %             /         \ /         \
    
19. %            A --------- C ----------E
    
20. %
    
21. %  problem3.m       Dianne P. O'Leary  12/2005
    
22. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
23. 
    
24. 
    
25. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
26. %
    
27. % The original problem  A f = b, where f is the force vector
    
28. % and b is the load vector.
    
29. % The variables (columns) are the unknown forces:
    
30. % The equations (rows) balance horizontal and vertical forces
    
31. %  at each node.
    
32. %
    
33. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
34. 
    
35. c = cos(pi/3);
    
36. s = sin(pi/3);
    
37. l = -50;
    
38. 
    
39. %     Ah Av AB BC CD DE BD AC CE Ev
    
40. A = [ -1  0  c  0  0  0  0  1  0  0   % horizontal forces at A
    
41.        0  1  s  0  0  0  0  0  0  0   % vertical   forces at A
    
42.        0  0 -c  c  0  0  1  0  0  0   % horizontal forces at B
    
43.        0  0  s  s  0  0  0  0  0  0   % vertical   forces at B
    
44.        0  0  0 -c  c  0  0 -1  1  0   % horizontal forces at C
    
45.        0  0  0  s  s  0  0  0  0  0   % vertical   forces at C
    
46.        0  0  0  0 -c  c -1  0  0  0   % horizontal forces at D
    
47.        0  0  0  0  s  s  0  0  0  0   % vertical   forces at D
    
48.        0  0  0  0  0 -c  0  0 -1  0   % horizontal forces at E
    
49.        0  0  0  0  0  s  0  0  0 -1]; % vertical   forces at E
    
50. 
    
51. b = [  0; 0; 0; 0; 0; l; 0; 0; 0; 0];
    
52. 
    
53. f = A \ b
    
54. 
    
55. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
56. %
    
57. %  Solve the modified problem: move node E so that member CE is
    
58. %  2 times the length of the other members.
    
59. %  The resulting matrix is A - Z V^T.
    
60. %
    
61. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
62. 
    
63. 
    
64. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
65. % Angle DCE remains 60 degrees.
    
66. % Angle CDE is 90 degrees
    
67. % Angle DEC is 30 degrees.
    
68. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
69. 
    
70. cosE =  cos(pi/6)
    
71. sinE =  sin(pi/6)
    
72. cosC =  cos(pi/3);
    
73. 
    
74. gamma = pi/6;   % = angle at D between member 4 and horizontal.
    
75. cosgam = cos(gamma);
    
76. singam = sin(gamma);
    
77. 
    
78. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
79. % New column 6 of matrix:
    
80. % Member 4 applies horizontal force at D with factor  cos gamma   
    
81. %                  vertical   force at D with factor  sin gamma   
    
82. % Member 4 applies horizontal force at E with factor -cos gamma   
    
83. %                  vertical   force at E with factor  sin gamma   
    
84. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
85. 
    
86. V = zeros(10,1);
    
87. V(6,1) = 1;
    
88. 
    
89. Z = [A(:,6) - [0 0 0 0   0      0   cosgam singam -cosgam singam]'];
    
90. 
    
91. [L,U,P]=lu(A);
    
92. 
    
93. fmod = sherman_mw(L,U,b,Z,V,P)
    
94. 
    
95. fmodtru = (A-Z*V')\b;
    
96. 
    
97. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
98. % Display the results.
    
99. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
100. 
     
101. disp(sprintf('The original vector of forces:'))
     
102. disp(f')
     
103. disp(sprintf('The vector of forces after the node is moved:'))
     
104. disp(fmod')
     
105. disp(sprintf('The difference between the Sherman-Morrison-Woodbury'))
     
106. disp(sprintf(' answer and the answer computed with backslash: %e',...
     
107.      norm(fmod-fmodtru)))

复制代码

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %  Solution to Problem 4 in the homework assignment,
   
3. %  Updating and Downdating Matrix Factorizations:
   
4. %     A Change in Plans
   
5. %
   
6. %  For matrices of various sizes n, find the number of updates
   
7. %  that makes the cost of using the Sherman-Morrison-Woodbury
   
8. %  formula to solve a modified linear system comparable to
   
9. %  solving the problem with backslash.
   
10. %
    
11. %  problem4.m       Dianne P. O'Leary  12/2005
    
12. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
13. 
    
14. nn = [10:10:50,100:100:1000];
    
15. 
    
16. ntries = 5;   % Number of repeats, used to try to make the timings
    
17.               % more reproducible.
    
18. k = 0;
    
19. 
    
20. for nsize=1:length(nn);
    
21. 
    
22. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
23. %  Generate a random problem of size n.
    
24. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
25. 
    
26.    n = nn(nsize);
    
27.    A = rand(n,n);
    
28.    x = rand(n,1);
    
29.    b = A*x;
    
30.    [L,U,P]=lu(A);
    
31.    P = sparse(P);
    
32. 
    
33. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
34. %  Set timefact to be the time for solving using backslash.
    
35. %  (Ignore the overhead in the "for" loop.)
    
36. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
37. 
    
38.    tic
    
39.    for j=1:ntries,
    
40.       xcomp = A \ b;
    
41.    end
    
42.    timefact(nsize) = toc/ntries;
    
43.    ttime = 0;
    
44. 
    
45. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
46. %  Use Sherman-Morrison-Woodbury with a rank k update, and increase
    
47. %  k until the time is greater than that for backslash.
    
48. %  (Ignore the overhead in the "for" loop.)
    
49. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
50. 
    
51.    while (ttime < timefact(nsize))&(k < n),
    
52.      k = k + 1;
    
53.      Z = rand(n,k);
    
54.      V = rand(n,k);
    
55. 
    
56.      tic
    
57.      for j=1:ntries
    
58.         xnew = sherman_mw(L,U,b,Z,V,P);
    
59.      end
    
60.      ttime = toc/ntries;
    
61.    end
    
62.    cross(nsize) = k;
    
63.    k = k - 1;
    
64. 
    
65. end
    
66. 
    
67. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
68. %  Plot the results.
    
69. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
70. 
    
71. plot(nn,cross,'LineWidth',1.5)
    
72. xlabel('Size of matrix')
    
73. ylabel('Rank of update')
    
74. title('Making Sherman-Morrison-Woodbury time comparable to Backslash')
    
75. print -dtiff smw.tif

复制代码

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %  Solution to Problem 6 in the homework assignment,
   
3. %  Updating and Downdating Matrix Factorizations:
   
4. %     A Change in Plans
   
5. %
   
6. %  We solve two truss problems, the first through matrix
   
7. %  factorization and the second using a QR update.
   
8. %
   
9. %  First truss:  All members are equal length.
   
10. %                Assume a load of 50 at C.
    
11. %
    
12. %                  B           D
    
13. %                 / \         / \
    
14. %                /   \       /   \
    
15. %               /     \     /     \
    
16. %              /       \   /       \
    
17. %             /         \ /         \
    
18. %            A --------- C ----------E
    
19. %
    
20. %
    
21. %  problem6.m       Dianne P. O'Leary  12/2005
    
22. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
23. 
    
24. 
    
25. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
26. %
    
27. %  The original problem  A f = b, where f is the force vector
    
28. %  and b is the load vector.
    
29. % The variables (columns) are the unknown forces:
    
30. % The equations (rows) balance horizontal and vertical forces
    
31. %  at each node.
    
32. %
    
33. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
34. 
    
35. c = cos(pi/3);
    
36. s = sin(pi/3);
    
37. l = -50;
    
38. 
    
39. %     Ah Av AB BC CD DE BD AC CE Ev
    
40. A = [ -1  0  c  0  0  0  0  1  0  0   % horizontal forces at A
    
41.        0  1  s  0  0  0  0  0  0  0   % vertical   forces at A
    
42.        0  0 -c  c  0  0  1  0  0  0   % horizontal forces at B
    
43.        0  0  s  s  0  0  0  0  0  0   % vertical   forces at B
    
44.        0  0  0 -c  c  0  0 -1  1  0   % horizontal forces at C
    
45.        0  0  0  s  s  0  0  0  0  0   % vertical   forces at C
    
46.        0  0  0  0 -c  c -1  0  0  0   % horizontal forces at D
    
47.        0  0  0  0  s  s  0  0  0  0   % vertical   forces at D
    
48.        0  0  0  0  0 -c  0  0 -1  0   % horizontal forces at E
    
49.        0  0  0  0  0  s  0  0  0 -1]; % vertical   forces at E
    
50. 
    
51. b = [  0; 0; 0; 0; 0; l; 0; 0; 0; 0];
    
52. 
    
53. f = A \ b
    
54. 
    
55. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
56. %
    
57. %  Solve the modified problem: move node E so that member CE is
    
58. %  2 times the length of the other members.
    
59. %  The resulting matrix is A - Z V^T.
    
60. %
    
61. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
62. 
    
63. 
    
64. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
65. % Angle DCE remains 60 degrees.
    
66. % Angle CDE is 90 degrees
    
67. % Angle DEC is 30 degrees.
    
68. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
69. 
    
70. cosE =  cos(pi/6)
    
71. sinE =  sin(pi/6)
    
72. cosC =  cos(pi/3);
    
73. 
    
74. gamma = pi/6;   % = angle at D between member 4 and horizontal.
    
75. cosgam = cos(gamma);
    
76. singam = sin(gamma);
    
77. 
    
78. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
79. % New column 6 of matrix:
    
80. % Member 4 applies horizontal force at D with factor  cos gamma   
    
81. %                  vertical   force at D with factor  sin gamma   
    
82. % Member 4 applies horizontal force at E with factor -cos gamma   
    
83. %                  vertical   force at E with factor  sin gamma   
    
84. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
85. 
    
86. V = zeros(10,1);
    
87. V(6,1) = 1;
    
88. 
    
89. Z = [A(:,6) - [0 0 0 0   0      0   cosgam singam -cosgam singam]'];
    
90. 
    
91. [L,U,P]=lu(A);
    
92. 
    
93. [Q,R]=qr(A);
    
94. 
    
95. global Anew
    
96. Anew = A - Z*V';
    
97. 
    
98. [Qnew,Rnew] = qrcolchange(Q,R,6,-Z)
    
99. 
    
100. fmod = Rnew\(Qnew'*b);
     
101. fmodtru = (A-Z*V')\b;
     
102. 
     
103. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
104. % Display the results.
     
105. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
106. 
     
107. disp(sprintf('The original vector of forces:'))
     
108. disp(f')
     
109. disp(sprintf('The vector of forces after the node is moved:'))
     
110. disp(fmod')
     
111. disp(sprintf('The difference between the QR-update'))
     
112. disp(sprintf(' answer and the answer computed with backslash: %e',...
     
113.      norm(fmod-fmodtru)))

复制代码

本帖隐藏的内容

here？Dianne P. O'Leary-Scientific Computing with Case Studies-Society for .pdf (4.52 MB)

2016-3-20 22:49:08 上传

我来吧  哦 已经有人应助了 没看见