Computational Science and Engineering课后部分编程习题解答（MATLAB实现）

1.3.9

题目

The Cholesky command A= choi (K) produces an upper triangular A with K = AT A. The square roots of the pivots from D are now included on the diagonal of A (so Cholesky fails unless K = KT and the pivots are positive) . Try the chol command on K3 , T3 , B3 , and B3 + eps * eye(3) .

分析

这题比较简单，只需要构建K3、T3、B3然后调用函数就可以了,从结果可以看出，没有执行chol(B3)，说明B3不是正定的，为了生成KTBC矩阵的方便，先写一个子函数用于生成KTBC四种矩阵，函数如下：

KTBC

function [K,T,B,C] = KTBC(n)
% Create the four special matrices assuming n>1
K = toeplitz ([2 -1 zeros(1,n-2)]);
T = K; T(1,1) = 1;
B = K; B(1,1) = 1; B(n,n) = 1;
C = K; C(1,n) = -1; C(n,1) = -1; 
end

hw_1_3_9

function hw_1_3_9()
[K3,T3,B3,~]=KTBC(3);
fprintf('the result of 1.3.9 is below\n');
fprintf('The result of chol(K3) is below:\n');
disp(chol(K3));
fprintf('****************************\n');
fprintf('The result of chol(T3) is below:\n');
disp(chol(T3));
fprintf('****************************\n');
try
fprintf('The result of chol(B3) is below:\n');
disp(chol(B3));
fprintf('****************************\n');
catch
fprintf('The result of chol(B3 + eps * eye(3)) is below:\n');
disp(chol(B3 + eps * eye(3)));
fprintf('****************************\n');
end
end

结果

the result of 1.3.9 is below
The result of chol(K3) is below:
    1.4142   -0.7071         0
         0    1.2247   -0.8165
         0         0    1.1547

****************************
The result of chol(T3) is below:
     1    -1     0
     0     1    -1
     0     0     1

****************************
The result of chol(B3) is below:
The result of chol(B3 + eps * eye(3)) is below:
    1.0000   -1.0000         0
         0    1.0000   -1.0000
         0         0    0.0000

****************************

1.5.4

题目

Continue 3 to find an eigenvector matrix Q by [Q, E] = eig(K) . The Discrete Sine Transform DST = Q * diag( [ -1 -1 1 -1 1 ] ) starts each column with a positive entry. The matrix J K = [1 : 5 ] ’ * [1 : 5] has entries j times k. Verify that DST agrees with sin( J K * pi/6) jsqrt(3) , and test DSTT = DST-1 .

1.5.33

题目

Scarymatlab shows what can happen when roundoff destroys symmetry: A = [ 1 1111; 1:5] ’ ; B=A’ * A; P=A * inv(B) * A’ ; [ Q, E] = eig§; B is exactly symmetric. The projection P should be symmetric, but isn’t. From Q’ * Q show that two eigenvectors of P fail badly to have inner product 0.

分析

这题有点奇怪，题目说P不是对称的，但是我的结果显示是对称的，可能是他写书的时候matlab的精度还没那么高，具体原因不太清楚。

hw_1_5_33

function hw_1_5_33()
A=[1 1 1 1 1;1:5]';
B=A'*A;
P=A*inv(B)*A';
[Q,E]=eig(P);
fprintf('the result of 1.3.9 is below\n');
fprintf('the result of P is below\n');
disp(P);
fprintf('the result of Q is below\n');
disp(Q);
fprintf('the result of E is below\n');
disp(E);
end

结果

the result of 1.3.9 is below
the result of P is below
    0.6000    0.4000    0.2000    0.0000   -0.2000
    0.4000    0.3000    0.2000    0.1000         0
    0.2000    0.2000    0.2000    0.2000    0.2000
    0.0000    0.1000    0.2000    0.3000    0.4000
   -0.2000         0    0.2000    0.4000    0.6000

the result of Q is below
    0.7746    0.6325   -0.0488   -0.2239   -0.1895
    0.5164   -0.6325    0.1497    0.7346    0.4605
    0.2582   -0.3162    0.3482   -0.5987    0.0665
    0.0000   -0.0000    0.5466   -0.1107   -0.7564
   -0.2582    0.3162    0.7451    0.1987    0.4189

the result of E is below
    1.0000         0         0         0         0
         0   -0.0000         0         0         0
         0         0    1.0000         0         0
         0         0         0   -0.0000         0
         0         0         0         0    0.0000

1.7.27

题目

For n = 3 and n = 4, find the QR factors of the special tridiagonal matrices T and K and B from qr(T) and qr(K) and qr(B) . Can you see a pattern?

分析

这题比较简单，但是用个循环可以比较简洁

hw_1_7_27

function hw_1_7_27()
names=['K','T','B'];
for i = 3:4
    [K,T,B,~]=KTBC(i);
    disp(juzhen(1));
    for j = 1:3
        [Q,R]=qr(juzhen(1:i,(j-1)*i+1:j*i));
        name=names(j)+string(i);
        fprintf('the Q of '+name+ ' is below\n');
        disp(Q);
        fprintf('the R of '+name+ ' is below\n');
        disp(R);
        fprintf('*************\n')
    end
end
end

结果

the Q of K3 is below
   -0.8944   -0.3586    0.2673
    0.4472   -0.7171    0.5345
         0    0.5976    0.8018

the R of K3 is below
   -2.2361    1.7889   -0.4472
         0   -1.6733    1.9124
         0         0    1.0690

*************
the Q of T3 is below
   -0.7071   -0.4082    0.5774
    0.7071   -0.4082    0.5774
         0    0.8165    0.5774

the R of T3 is below
   -1.4142    2.1213   -0.7071
         0   -1.2247    2.0412
         0         0    0.5774

*************
the Q of B3 is below
   -0.7071   -0.4082    0.5774
    0.7071   -0.4082    0.5774
         0    0.8165    0.5774

the R of B3 is below
   -1.4142    2.1213   -0.7071
         0   -1.2247    1.2247
         0         0   -0.0000

*************
the Q of K4 is below
   -0.8944   -0.3586   -0.1952    0.1826
    0.4472   -0.7171   -0.3904    0.3651
         0    0.5976   -0.5855    0.5477
         0         0    0.6831    0.7303

the R of K4 is below
   -2.2361    1.7889   -0.4472         0
         0   -1.6733    1.9124   -0.5976
         0         0   -1.4639    1.9518
         0         0         0    0.9129

*************
the Q of T4 is below
   -0.7071   -0.4082   -0.2887    0.5000
    0.7071   -0.4082   -0.2887    0.5000
         0    0.8165   -0.2887    0.5000
         0         0    0.8660    0.5000

the R of T4 is below
   -1.4142    2.1213   -0.7071         0
         0   -1.2247    2.0412   -0.8165
         0         0   -1.1547    2.0207
         0         0         0    0.5000

*************
the Q of B4 is below
   -0.7071   -0.4082   -0.2887    0.5000
    0.7071   -0.4082   -0.2887    0.5000
         0    0.8165   -0.2887    0.5000
         0         0    0.8660    0.5000

the R of B4 is below
   -1.4142    2.1213   -0.7071         0
         0   -1.2247    2.0412   -0.8165
         0         0   -1.1547    1.1547
         0         0         0   -0.0000

*************

2.3.3

题目

A QR experiment on the Vandermonde matrix V = fl iplr(vander((O: 49)/49) ) is proposed by Trefethen-Bau [159] . A= V(:, 1: 12) and b = cos(O: .08: 3.92) 1 have m = 50 and n = 12. In format long compute u in multiple ways (MATLAB’s svd(A) is another way) . How many correct digits does each method yield?

1. Directly from the normal equations by AT A(ATb)
2. From the unmodified Gram-Schmidt codes by R(QTb)
3. From the modified Gram-Schmidt code with v in line 4
4. 4. From the Householder code using 12 columns of Q and 12 rows of R 5. From MATLAB’s A\b and MATLAB’s qr (which uses Householder)

3.2.19

题目

就是拟合一下hat function的曲线（复制过来有点bug就直接中文描述了）

分析

比较简单，就是spline函数的使用，可以看出点越多越接近真实的hat函数。

hat

function y = hat(x)
y=x.*(x<0.5)+(1-x).*(x>=0.5);
end

hw_3_2_19

function hw_3_2_19()
for i=1:2
z=3^i;
x=0:1/z:1;
y=hat(x);
xx=0:.025:1;
yy=spline(x,y,xx);
plot(x,y,'o',xx,yy);
hold on;
end
end

结果

3.2.20

题目

和上一题差不多，把hat换成step而已

mystep

function y = mystep(x)
y=0.*(x<0.5)+1.*(x>=0.5);
end

hw_3_2_20

function hw_3_2_20()
for i=1:2
z=3^i;
x=0:1/z:1;
y=mystep(x);
xx=0:.025:1;
yy=spline(x,y,xx);
plot(x,y,'o',xx,yy);
hold on;
end
end

结果

p.s. 有新的会更