Chapter 6. MATLAB数值计算基础

课后习题解答

转载请注明出处： http://blog.csdn.net/i_love_home/article/details/17514969

还在学习中，一起交流

1. 利用MATLAB提供的randn函数生成符合正态分布的10*5随机矩阵A，进行如下操作。

（1）A各列元素的均值和标准方差。

（2）A的最大元素和最小元素。

（3）求A每行元素的和以及全部元素之和。

（4）分别对A的每列元素按升序、每行元素按降序排序。

A = randn(10, 5);
  % (1)
mean(A(:)); 
var(A(:));
  % (2)
max(A(:));
min(A(:));
  % (3)
sum(A, 2);
sum(A(:));
  % (4)
[Y1, I1] = sort(A, 1, 'ascend');
[Y2, I2] = sort(A, 2, 'descend');

2. 按要求对指定函数进行插值和拟合。

（1）按表6.4用3次样条方法插值计算0~90度范围内整数点的正弦值和0~75度范围内整数点的正切值，然后用5次多项式拟合方法计算相同的函数值，并将两种计算结果进行比较。

arf = (0 : 15 : 75)/180*pi; 
sarf = [0, 0.2588, 0.5000, 0.7071, 0.8660, 0.9659, 1.0000];
sarf(end) = [];
tarf = [0, 0.2679, 0.5774, 1.0000, 1.7320, 3.7320];
t = (0:0.1:75)/180*pi;
cs3 = interp1(arf, sarf, t, 'spline');
ct3 = interp1(arf, tarf, t, 'spline');
ps = polyfit(arf, sarf, 5);
pt = polyfit(arf, tarf, 5);
cs5 = polyval(ps, t);
ct5 = polyval(pt, t);
plot(t, cs3, 'b:', t, cs5, 'r-', t, ct3, 'b:', t, ct5, 'r-');
grid on; axis equal;
xlabel('x-axis'); ylabel('y-axis'); title('2-1');

（2）按表6.5用3次多项式方法插值计算1~100之间整数的平方根。

n = 1 : 10;
N = n.^2;
t = 1 : 0.4 : 10;
p = polyfit(n, N, 3);
f = polyval(p, t);
plot(n, N, ':o', t, f, '-*');
grid on; axis equal;
xlabel('x-axis'); ylabel('y-axis'); title('2-2');

3. 已知一组实验数据如表6.6所示。

n = 1 : 5;
x = [165, 123, 150, 123, 141];
y = [187, 126, 172, 125, 148];
[x, I] = sort(x, 2, 'ascend');
ty = [];
for i = 1 : size(I, 2)
    ty = [ty, y(I(i))];
end
y = ty;
X = min(x) : 0.5 : max(x); 
P = polyfit(x, y, 2);
Y = polyval(P, X);
plot(x, y, 'b', X, Y, 'r');
grid on; axis equal;
xlabel('x-axis'); ylabel('y-axis'); title('3');

4. 已知多项式

P1 = [0, 3, 2];
P2 = [5, -1, 2];
P3 = [1, 0, -0.5];
  % (1)
P = P1 + P2 + P3;
  % (2)
x = roots(P);
% poly(x)*6
  % (3)
x = 0.2 * (0:10);
y = polyval(P, x);

5. 求函数在指定区间的极值。

[U, fmin] = fminsearch('pro5_1', [1, 1]);

function f = pro5_1( u )
x = u(1); 
y = u(2);
f = 3*x.^2 + 2*x.*y + y.^2;

[x, fmin] = fminbnd('pro5_2', 0, pi);
fmin = -fmin;

function f = pro5_2( x )
f = sin(x) + cos(x^2);
f = -f;

6. 求函数在指定点的数值导数。

f = inline('sin(x).^2 + cos(x).^2');
x = 0 : 0.01 : 2*pi;
p = polyfit(x, f(x), 5);
dp = polyder(p);
dpx = polyval(dp, [pi/6, pi/4, pi/3, pi/2]);

  % (2) three solution
f = inline('sqrt(x.^2+1)');
g = inline('x./sqrt(x.^2+1)');
x = 0 : 0.01 : 4;
 
p = polyfit(x, f(x), 5);
dp = polyder(p);
dpx = polyval(dp, [1, 2, 3])
 
dx = diff(f([x, 4.01]))/0.01;
dg = g(x);
plot(x, dx, [1, 2, 3], dpx, 'o', x, dg, 'r');

7. 求数值积分

 % (1)
x = linspace(0, pi, 50);
y = sin(x).^5.*sin(5*x);
I = trapz(x, y);
 
f = inline('sin(x).^5.*sin(5*x)');
I = quadl(f, 0, pi);
 
I = quadl('pro7_1', 0, pi);

function f = pro7_1( x )
f = sin(x).^5.*sin(5*x);

  % (2)
f = inline('(1+x.^2)./(1+x.^4)');
I = quadl(f, -1, 1);

  % (3)
f = inline('x.*sin(x)./(1+cos(x).^2)');
quadl(f, 0, pi);

 % (4)
I = dblquad('pro7_4', 0, pi, 0, pi)

function f = pro7_4( x, y)
f = abs(cos(x+y));

f = inline('abs(cos(x+y))', 'x', 'y');
I = dblquad(f, 0, pi, 0, pi)

8. 已知h(t)=e^(-t),t≥0，取 N=64，对t从0~5s采样，用FFT做快速傅里叶变换，并绘制相应的振幅-频率图。

N = 64;
T = 5;
t = linspace(0, T, N);
h = exp(-t);
dt = t(2) - t(1);
f = 1/dt;
H = fft(h);
F = H(1:N/2+1);
f = f*(0:N/2)/N;
plot(f, abs(F), '-p');
grid on; axis equal; 
xlabel('Frequency'); ylabel('|F(k)|');

9. 分别用矩阵求逆、矩阵除法以及矩阵分解求线性方程组的解。

  % (1)
A = [2, 3, 5; 
     3, 7, 4; 
     1, -7, 1];
b = [10, 3, 5]';
 
X = A\b;
 
X = inv(A)*b;
 
[L, U] = lu(A);
X = U\(L\b);

  % (2)
A = [5, 1, -1, 0; 
     1, 0, 3, -1;
     -1, -1, 0, 5; 
     0, 0, 2, 4];
b = [1, 2, 3, -1]';
 
X = A\b;
 
X = inv(A)*b;
 
[L, U] = lu(A);
X = U\(L\b);

10. 求非齐次线性方程组的通解。

format rat;
A = [2, 1, -1, 1;
     4, 2, -2, 1;
     2, 1, -1, -1];
b = [1, 2, 1]';
[x, y] = LEGS(A, b);
format short;

function [x, y] = LEGS( A, b )
[row, col] = size(A);
y = [];
if norm(b) > 0
    if rank(A) == rank([A, b])
        if rank(A) == 0
            disp('原方程组有唯一解 x: ');
            x = A\b;
        else
            disp('原方程组有无穷个解，特解 x，其齐次方程组的基础解系 y');
            x = A\b;
            y = null(A, 'r');
        end
    else
        disp('方程组无解');
        x = [];
    end
else
    disp('原方程组有零解 x：');
    x = zeros(col, 1);
    if rank(A) < col
        diap('方程组有无穷个解，基础解系为 y');
        y = null(A, 'r');
    end
end

11. 求下列方程的根。

  % (1)
f = inline('x-sin(x)./x');
res = fzero(f, 0.5);
 
res = fzero('pro11_1', 0.5);

function f = pro11_1( x )
f = x - sin(x)/x;

  % (2)
f = inline('sin(x).^2.*exp(-0.1*x)-0.5*abs(x)');
res = fzero(f, 1.5);

12. 求非线性组在 (0.5, 0.5) 附近的数值解。

x = fsolve('pro12', [0.5, 0.5], optimset('Display', 'off'));

function f = pro12( par )
x = par(1);
y = par(2);
f(1) = x - 0.6*sin(x) - 0.3*cos(y);
f(2) = y - 0.6*cos(x) + 0.3*sin(y);

15. 已知线性方程组 ，其中

运用稀疏存储矩阵的方式求其解。

A = [2, -1, 0, 0, 0;
     -1, 2, -1, 0, 0;
     0, -1, 2, -1, 0;
     0, 0, -1, 2, -1;
     0, 0, 0, -1, 2];
A = sparse(A);
b = [1, 0, 0, 0, 0]';
X = A\b;