数值分析作业（第二章）：代码+手写计算

MATLAB_0">《数值计算方法》丁丽娟-数值实验作业-第二章（MATLAB）

作业P58: 1 ，2，3，6，8(1), 12, 13
数值实验P61: 2, 3

数值实验（第二章）

代码仓库：https://github.com/sylvanding/bit-numerical-analysis-hw

P61. 2

试用MATLAB软件编程实现追赶法求解三对角方程组的算法，并考虑如下梯形电阻电路问题：

其中电路中的各个电流 $\{i_1, i_2, \dots, i_8\}$ 满足下列线性方程组：
$\begin{aligned} 2 i_1-2 i_2 &=V / R \\ -2 i_1+5 i_2-2 i_3&=0 \\ -2 i_2+5 i_3-2 i_4&=0 \\ -2 i_3+5 i_4-2 i_5&=0 \\ -2 i_4+5 i_5-2 i_6&=0 \\ -2 i_5+5 i_6-2 i_7&=0 \\ -2 i_6+5 i_7-2 i_8&=0 \\ -2 i_7+5 i_8&=0 \end{aligned}$
设 $V=220\text{V}, R=27\Omega$ ，求各段电路的电流量。

实验内容、步骤及结果

chap-2\thomas_algorithm.m:

function x = thomas_algorithm(a, b, c, d)% a: sub-diagonal (length n-1)% b: main diagonal (length n)% c: super-diagonal (length n-1)% d: right-hand side (length n)n = length(b);% Forward eliminationfor i = 2:nw = a(i-1) / b(i-1);b(i) = b(i) - w * c(i-1);d(i) = d(i) - w * d(i-1);end% Back substitutionx = zeros(n, 1);x(n) = d(n) / b(n);for i = n-1:-1:1x(i) = (d(i) - c(i) * x(i+1)) / b(i);end
end

chap-2\Q2.m:

% Define input parameters
a = -2;
a = repmat(a, 1, 7);
b = [2, 5, 5, 5, 5, 5, 5, 5];
c = a;
V = 220;
R = 27;
d = [V/R, 0, 0, 0, 0, 0, 0, 0];% Call the function
x = thomas_algorithm(a, b, c, d);% Display the output
disp(x);

实验结果分析

>> Q28.14784.07372.03651.01750.50730.25060.11940.0477

所求电阻输出如上。

P61. 3

方程组的性态和矩阵的条件数的实验。设有线性方程组 $A x = b$ ，其中系数矩阵 $A=(a_{ij})_{n\times n}$ 的元素分别为 $a_{i j}=x_i^{j-1}\left(x_i=1+0.1 i, i, j=1,2, \cdots, n\right)$ 或 $a_{i j}=\frac{1}{i+j-1}(i, j=1,2, \cdots, n)$ ，右端向量 $\boldsymbol{b}=\left(\sum_{j=1}^n a_{1 j}, \sum_{j=1}^n a_{2 j}, \cdots, \sum_{j=1}^n a_{n j}\right)^{\mathrm{T}}$ . 利用MATLAB中的库函数，

取 $n = 5, 10, 20$ ，分别求出系数矩阵的2-条件数，判别它们是否是病态阵？随着 $n$ 的增大，矩阵性态的变化如何？
分别取 $n = 5, 10, 20$ ，解上述两个线性方程组 $A x = b$ ，并将求得的解与精确解 $\cdots, 1)^{\mathrm{T}}$ 作比较，说明了什么？
取 $n = 10$ ，对系数矩阵中的 $a_{22}$ 和 $a_{nn}$ 增加一个扰动 $10^{-8}$ ，求解方程组 $(A+\delta A) x=b$ ，解的变化如何？

实验内容、步骤及结果

chap-2\generateLinearSystem.m:

function [A, b] = generateLinearSystem(n, matrixType)% This function generates the coefficient matrix A and the vector b for% the linear system Ax = b. The parameter n specifies the size of the% matrix, and matrixType specifies the type of matrix to generate.%% matrixType can be:% 'polynomial' - A is generated using a_{i j} = x_i^{j-1} where x_i = 1 + 0.1 * i% 'hilbert'    - A is generated using a_{i j} = 1 / (i + j - 1)if strcmp(matrixType, 'polynomial')% Generate polynomial matrixA = zeros(n);for i = 1:nx_i = 1 + 0.1 * i;for j = 1:nA(i, j) = x_i^(j-1);endendelseif strcmp(matrixType, 'hilbert')% Generate Hilbert matrixA = zeros(n);for i = 1:nfor j = 1:nA(i, j) = 1 / (i + j - 1);endendelseerror('Invalid matrix type. Use ''polynomial'' or ''hilbert''.');end% Compute vector bb = sum(A, 2);
end

chap-2\Q31.m:

% Define the values of n
n_values = [5, 10, 20];
matrixType = 'polynomial'; % or 'hilbert'% Loop through each value of n
for n = n_values% Generate the linear system[A, b] = generateLinearSystem(n, matrixType);% Compute the 2-norm condition numbercond_num = cond(A, 2);% Determine if the matrix is ill-conditionedif cond_num > 1e10 % Cond(A) >> 1result = 'Ill-conditioned';elseresult = 'Well-conditioned';end% Display the condition number and the resultfprintf('%s, n = %d, 2-norm condition number = %.2e, Result: %s\n', matrixType, n, cond_num, result);
end

chap-2\Q32.m:

% Define the values of n
n_values = [5, 10, 20];
matrixTypes = {'polynomial', 'hilbert'};% Exact solution
x_exact = @(n) ones(n, 1);% Loop through each matrix type and value of n
for matrixType = matrixTypesfor n = n_values% Generate the linear system[A, b] = generateLinearSystem(n, matrixType{1});% Solve the linear systemx_approx = A \ b;% Compute the relative error using the infinity normrel_error = norm(x_approx - x_exact(n), inf) / norm(x_exact(n), inf);% Display the relative errorfprintf('%s, n = %d, Relative Error (inf-norm) = %.2e\n', matrixType{1}, n, rel_error);end
end

chap-2\Q33.m:

% Define the value of n
n = 10;
matrixType = 'polynomial'; % or 'hilbert'% Generate the original linear system
[A, b] = generateLinearSystem(n, matrixType);% Exact solution
x_exact = @(n) ones(n, 1);% Solve the original system
x_original = A \ b;% Introduce perturbation
A_perturbed = A;
A_perturbed(2, 2) = A_perturbed(2, 2) + 1e-8;
A_perturbed(n, n) = A_perturbed(n, n) + 1e-8;% Solve the perturbed system
x_perturbed = A_perturbed \ b;% Compute relative errors
relative_error_exact = norm(x_original - x_exact(n), inf) / norm(x_exact(n), inf);
relative_error_perturbed = norm(x_perturbed - x_exact(n), inf) / norm(x_exact(n), inf);% Display the results
fprintf('%s\n', matrixType);
fprintf('Relative error of the original solution: %.2e\n', relative_error_exact);
fprintf('Relative error of the perturbed solution: %.2e\n', relative_error_perturbed);

实验结果分析

Q31

>> Q31
polynomial, n = 5, 2-norm condition number = 5.36e+05, Result: Well-conditioned
polynomial, n = 10, 2-norm condition number = 8.68e+11, Result: Ill-conditioned
polynomial, n = 20, 2-norm condition number = 6.07e+21, Result: Ill-conditioned
>> Q31
hilbert, n = 5, 2-norm condition number = 4.77e+05, Result: Well-conditioned
hilbert, n = 10, 2-norm condition number = 1.60e+13, Result: Ill-conditioned
hilbert, n = 20, 2-norm condition number = 2.11e+18, Result: Ill-conditioned

随 $n$ 增大，系数矩阵的病态程度进一步提升。

Q32

>> Q32
polynomial, n = 5, Relative Error (inf-norm) = 2.89e-11
polynomial, n = 10, Relative Error (inf-norm) = 7.21e-05
polynomial, n = 20, Relative Error (inf-norm) = 3.50e+05
hilbert, n = 5, Relative Error (inf-norm) = 6.17e-13
hilbert, n = 10, Relative Error (inf-norm) = 3.18e-04
hilbert, n = 20, Relative Error (inf-norm) = 2.64e+01Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  3.713214e-23.
> In Q32 (line 15): x_approx = A \ b;
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.276108e-19.
> In Q32 (line 15): x_approx = A \ b;

随 $n$ 增大，系数矩阵的病态程度进一步提升，造成解的相对误差进一步增大，且这个误差与方程组的性态有关。

Q33

>> Q33
polynomial
Relative error of the original solution: 7.21e-05
Relative error of the perturbed solution: 3.16e-01
>> Q33
hilbert
Relative error of the original solution: 3.18e-04
Relative error of the perturbed solution: 9.77e+00

当 $n = 10$ 时，线性系统polynomial和hilbert均为病态，即使引入一个为小扰动，也会产生较大的相对误差。