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

news/2024/10/3 18:12:14 标签: 数值计算, MATLAB, 北京理工大学

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 , … , 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 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

方程组的性态和矩阵的条件数的实验。设有线性方程组 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=1n​a1j​,∑j=1n​a2j​,⋯,∑j=1n​anj​)T. 利用MATLAB中的库函数，

  1. 取 n = 5 , 10 , 20 n=5,10,20 n=5,10,20，分别求出系数矩阵的2-条件数，判别它们是否是病态阵？随着 n n n的增大，矩阵性态的变化如何？
  2. 分别取 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作比较，说明了什么？
  3. 取 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 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 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+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 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均为病态，即使引入一个为小扰动，也会产生较大的相对误差。

手写作业

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述


http://www.niftyadmin.cn/n/5688788.html

相关文章

使用 Vue3 和 Axios 实现 CRUD 操作

使用 Vue3 和 Axios 实现 CRUD 操作

文章目录 1、准备工作2、创建 Vue 3 项目3、项目结构4、实现 CRUD 操作5、运行项目6、小结在当今的前端开发中，Vue.js 作为一款流行的 JavaScript 框架，正在被越来越多的开发者所青睐。尤其是 Vue 3 引入了 Composition API 和更优雅的响应式处理，使得模板编写和状态管理变得…
阅读更多...
Python、C++、java阶乘算法

Python、C++、java阶乘算法

最近，我除了Python还学了C和Java，然后在网上看到编程考题：阶乘。 首先，我们先理解什么是阶乘。 阶乘是数学中的一个概念，通常定义为从1乘到指定的数。具体来说，一个正整数的阶乘（记作n!&#…
阅读更多...
RAC被修改权限及相关问题

RAC被修改权限及相关问题

RDBMS ： 19.19 修改RAC权限及相关问题 修改RAC权限，参考文档： How to check and fix file permissions on Grid Infrastructure environment (Doc ID 1931142.1) Script to capture and restore file permission in a directory (for eg. O…
阅读更多...
接口 抽象类

接口 抽象类

接口和抽象类都是用来实现面向对象编程中的抽象概念的工具。 接口是一种抽象的数据类型，它定义了一组抽象方法。接口中的方法没有具体的实现，只有方法的声明。类可以实现一个或多个接口，并实现接口中的方法。接口提供了一种规范，…
阅读更多...
uniapp实战教程：如何封装一个可复用的表单组件

uniapp实战教程：如何封装一个可复用的表单组件

在uniapp开发过程中，表单组件的使用场景非常广泛。为了提高开发效率，我们可以将常用的表单组件进行封装。本文将带你了解如何在uniapp中封装一个表单组件，让你只需要通过属性配置轻松实现各种表单，效果图如下： 一、准备…
阅读更多...
JavaScript break与continue语句

JavaScript break与continue语句

break语句和continue语句都具有跳转作用&#xff0c;可以让代码不按既有的顺序执行。 break break语句用于跳出代码块或循环 for(i0;i<100;i){if(i5){break;}console.log(i);} continue continue语句用于应即终止本轮循环,返回循环结构的头部&#xff0c;开始下一轮循环。…
阅读更多...
【Kubernetes】常见面试题汇总（五十四）

【Kubernetes】常见面试题汇总（五十四）

目录 120.创建 init C 容器后&#xff0c;其状态不正常&#xff1f; 特别说明&#xff1a; 题目 1-68 属于【Kubernetes】的常规概念题&#xff0c;即 “ 汇总&#xff08;一&#xff09;~&#xff08;二十二&#xff09;” 。 题目 69-113 属于【Kubernetes】的生产…
阅读更多...
Hive数仓操作（八）

Hive数仓操作（八）

一、Hive中的分桶表 1. 分桶表的概念 分桶表是Hive中一种用于提升查询效率的表类型。分桶指的是根据指定列的哈希值将数据划分到不同的文件&#xff08;桶&#xff09;中。 2. 分桶表的原理 哈希分桶&#xff1a;根据分桶列计算哈希值&#xff0c;对哈希值取模&#xff0c;将…
阅读更多...
最新文章

Copyright @ 2022~2023