《数值计算方法》丁丽娟-数值实验作业-第二章(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 , … , i 8 } \{i_1, i_2, \dots, i_8\} {i1,i2,…,i8}满足下列线性方程组:
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 \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} 2i1−2i2−2i1+5i2−2i3−2i2+5i3−2i4−2i3+5i4−2i5−2i4+5i5−2i6−2i5+5i6−2i7−2i6+5i7−2i8−2i7+5i8=V/R=0=0=0=0=0=0=0
设 V = 220 V , R = 27 Ω V=220\text{V}, R=27\Omega V=220V,R=27Ω,求各段电路的电流量。
实验内容、步骤及结果
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 Ax=b Ax=b,其中系数矩阵 A = ( a i j ) n × n A=(a_{ij})_{n\times n} A=(aij)n×n的元素分别为 a i j = x i j − 1 ( x i = 1 + 0.1 i , i , j = 1 , 2 , ⋯ , n ) a_{i j}=x_i^{j-1}\left(x_i=1+0.1 i, i, j=1,2, \cdots, n\right) aij=xij−1(xi=1+0.1i,i,j=1,2,⋯,n)或 a i j = 1 i + j − 1 ( i , j = 1 , 2 , ⋯ , n ) a_{i j}=\frac{1}{i+j-1}(i, j=1,2, \cdots, n) aij=i+j−11(i,j=1,2,⋯,n),右端向量 b = ( ∑ j = 1 n a 1 j , ∑ j = 1 n a 2 j , ⋯ , ∑ j = 1 n a n j ) T \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}} b=(∑j=1na1j,∑j=1na2j,⋯,∑j=1nanj)T. 利用MATLAB中的库函数,
- 取 n = 5 , 10 , 20 n=5,10,20 n=5,10,20,分别求出系数矩阵的2-条件数,判别它们是否是病态阵?随着 n n n的增大,矩阵性态的变化如何?
- 分别取 n = 5 , 10 , 20 n=5,10,20 n=5,10,20,解上述两个线性方程组 A x = b Ax=b Ax=b,并将求得的解与精确解 x = ( 1 , 1 , ⋯ , 1 ) T x=(1,1, \cdots, 1)^{\mathrm{T}} x=(1,1,⋯,1)T作比较,说明了什么?
- 取 n = 10 n=10 n=10,对系数矩阵中的 a 22 a_{22} a22和 a n n a_{nn} ann增加一个扰动 1 0 − 8 10^{-8} 10−8,求解方程组 ( A + δ A ) x = b (A+\delta A) x=b (A+δ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 n 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 n 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 n=10 n=10时,线性系统polynomial和hilbert均为病态,即使引入一个为小扰动,也会产生较大的相对误差。