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

《数值计算方法》丁丽娟-数值实验作业-第二章（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​}满足下列线性方程组：
$$\begin{array}{rl}2i_1-2i_2& =V/R\\ -2i_1+5i_2-2i_3& =0\\ -2i_2+5i_3-2i_4& =0\\ -2i_3+5i_4-2i_5& =0\\ -2i_4+5i_5-2i_6& =0\\ -2i_5+5i_6-2i_7& =0\\ -2i_6+5i_7-2i_8& =0\\ -2i_7+5i_8& =0\end{array}$$
设$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 elimination
    for i = 2:n
        w = a(i-1) / b(i-1);
        b(i) = b(i) - w * c(i-1);
        d(i) = d(i) - w * d(i-1);
    end
    
    % Back substitution
    x = zeros(n, 1);
    x(n) = d(n) / b(n);
    for i = n-1:-1:1
        x(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);

实验结果分析

>> Q2
    8.1478
    4.0737
    2.0365
    1.0175
    0.5073
    0.2506
    0.1194
    0.0477

所求电阻输出如上。

P61. 3

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

1. 取$n=5,10,20$，分别求出系数矩阵的2-条件数，判别它们是否是病态阵？随着$n$的增大，矩阵性态的变化如何？
2. 分别取$n=5,10,20$，解上述两个线性方程组$Ax=b$，并将求得的解与精确解$x=(1,1,\cdots ,1{)}^{\mathrm{T}}$作比较，说明了什么？
3. 取$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 matrix
        A = zeros(n);
        for i = 1:n
            x_i = 1 + 0.1 * i;
            for j = 1:n
                A(i, j) = x_i^(j-1);
            end
        end
    elseif strcmp(matrixType, 'hilbert')
        % Generate Hilbert matrix
        A = zeros(n);
        for i = 1:n
            for j = 1:n
                A(i, j) = 1 / (i + j - 1);
            end
        end
    else
        error('Invalid matrix type. Use ''polynomial'' or ''hilbert''.');
    end

    % Compute vector b
    b = 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 number
    cond_num = cond(A, 2);
    
    % Determine if the matrix is ill-conditioned
    if cond_num > 1e10 % Cond(A) >> 1
        result = 'Ill-conditioned';
    else
        result = 'Well-conditioned';
    end
    
    % Display the condition number and the result
    fprintf('%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 = matrixTypes
    for n = n_values
        % Generate the linear system
        [A, b] = generateLinearSystem(n, matrixType{1});
        
        % Solve the linear system
        x_approx = A \ b;
        
        % Compute the relative error using the infinity norm
        rel_error = norm(x_approx - x_exact(n), inf) / norm(x_exact(n), inf);
        
        % Display the relative error
        fprintf('%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+01

Warning: 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均为病态，即使引入一个为小扰动，也会产生较大的相对误差。

手写作业

原文地址：https://blog.csdn.net/IYXUAN/article/details/142683390

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！