matlab使用教程(6)—线性方程组的求解

news2024/9/22 11:31:20
        进行科学计算时,最重要的一个问题是对联立线性方程组求解。在矩阵表示法中,常见问题采用以下形式:给定两个矩阵 A 和 b,是否存在一个唯一矩阵 x 使 Ax = b 或 xA = b?
        考虑一维示例具有指导意义。例如,方程
        7x = 21
        是否具有唯一解?
        答案当然是肯定的。方程有唯一解 x = 3。通过除法很容易求得该解:
        x = 21/7 = 3。
        该解通常不是通过计算 7 的倒数求得的,即先计算 7 –1 = 0.142857...,然后将 7 –1 乘以 21。这将需要更多的工作,而且如果 7 –1 以有限位数表示时,准确性会较低。类似注意事项也适用于多个未知数的线性方程组;MATLAB 在解此类方程时不会计算矩阵的逆。
        尽管这不是标准的数学表示法,但 MATLAB 使用标量示例中常见的除法术语来描述常规联立方程组的解。斜杠 / 和反斜杠 \ 这两个除号分别对应 MATLAB 函数 mrdivide mldivide。两种运算符分别用于未知矩阵出现在系数矩阵左侧或右侧的情况:
x = b/A
        表示使用 mrdivide 获得的矩阵方程 xA = b 的解。
x = A\b
        表示使用 mldivide 获得的矩阵方程 Ax = b 的解。
        考虑将方程 Ax = b 或 xA = b 的两端“除以”A。系数矩阵 A 始终位于“分母”中。
        x = A\b 的维度兼容性条件要求两个矩阵 A b 的行数相同。这样,解 x 的列数便与 b 的列数相同,并且其行维度等于 A 的列维度。对于 x = b/A ,行和列的角色将会互换。
        实际上,Ax=b 形式的线性方程组比 xA=b 形式的线性方程组更常见。因此,反斜杠的使用频率要远高于斜杠的使用频率。本节其余部分将重点介绍反斜杠运算符;斜杠运算符的对应属性可以从以下恒等式推知:
(b/A)' = (A'\b').
        系数矩阵 A 不需要是方阵。如果 A 的大小为 m×n,则有三种情况:
        m = n
        方阵方程组。求精确解。
        m > n
        超定方程组,即方程个数多于未知数个数。求最小二乘解。
        m < n
        欠定方程组,即方程个数少于未知数个数。使用最多 m 个非零分量求基本解。

1mldivide 算法

        mldivide 运算符使用不同的求解器来处理不同类型的系数矩阵。通过检查系数矩阵自动诊断各种情况。
        线性方程组 Ax = b 的通解描述了所有可能的解。可以通过以下方法求通解:
        1求对应的齐次方程组 Ax = 0 的解。使用 null 命令通过键入 null(A) 来执行此操作。这会将解空间的基向量恢复为 Ax = 0。任何解都是基向量的线性组合。
        2求非齐次方程组 Ax = b 的特定解。然后,可将 Ax = b 的任何解写成第 2 步中的 Ax = b 的特定解加上第 1 步中的基向量的线性组合之和。本节其余部分将介绍如何使用 MATLAB 求 Ax = b 的特定解,如第 2 步中所述。

2方阵方程组

        最常见的情况涉及到一个方阵系数矩阵 A 和一个右侧单列向量 b

2.1非奇异系数矩阵

        如果矩阵 A 是非奇异矩阵,则解 x = A\b 的大小与 b 的大小相同。例如:
A = pascal(3);
u = [3; 1; 4];
x = A\u
x =
10
-12
5
        可以确认 A*x 恰好等于 u。如果 A b 为方阵并且大小相同,则 x= A\b 也具有相同大小:
b = magic(3);
X = A\b
X =
19 -3 -1
-17 4 13
6 0 -6
        可以确认 A*x 恰好等于 b。以上两个示例具有确切的整数解。这是因为系数矩阵选为 pascal(3) ,这是满秩矩阵(非奇异的)。

2.2奇异系数矩阵

        如果方阵 A 不包含线性无关的列,则该矩阵为奇异矩阵。如果 A 为奇异矩阵,则 Ax = b 的解将不存在或不唯一。如果 A 接近奇异或检测到完全奇异性,则反斜杠运算符 A\b 会发出警告。如果 A 为奇异矩阵并且 Ax = b 具有解,可以通过键入以下内容求不是唯一的特定解
P = pinv(A)*b
        pinv(A) 是 A 的伪逆。如果 Ax = b 没有精确解,则 pinv(A) 将返回最小二乘解。例如:
A = [ 1 3 7 ; -1 4 4 ; 1 10 18 ] 为奇异矩阵,可以通过键入以下内容进行验证:
rank(A)
ans =
2
        由于 A 不是满秩,它有一些等于零的奇异值。对于 b =[5;2;12] ,方程 Ax = b 具有精确解:
pinv(A)*b
ans =
0.3850
-0.1103
0.7066
        通过键入以下内容验证 pinv(A)*b 是否为精确解
A*pinv(A)*b
ans =
5.0000
2.0000
12.0000
        但是,如果 b = [3;6;0] ,则 Ax = b 没有精确解。在这种情况下, pinv(A)*b 会返回最小二
乘解。键入
A*pinv(A)*b
ans =
-1.0000
4.0000
2.0000
        则不会返回原始向量 b。通过得到增广矩阵 [A b] 的简化行阶梯形式,可以确定 Ax = b 是否具有精确解。为此,对于此示例请输入
rref([A b])
ans =
1.0000 0 2.2857 0
0 1.0000 1.5714 0
0 0 0 1.0000
        由于最下面一行全部为零(最后一项除外),因此该方程无解。在这种情况下,pinv(A) 会返回最小二乘解。

3超定方程组

        此示例说明在对试验数据的各种曲线拟合中通常会如何遇到超定方程组。在多个不同的时间值 t 对数量 y 进行测量以生成以下观测值。可以使用以下语句输入数据并在表中查看该数据。
t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)
B=6×2 table
t y
___ ____
0 0.82
0.3 0.72
0.8 0.63
1.1 0.6
1.6 0.55
2.3 0.5
        尝试使用指数衰减函数对数据进行建模
 
        上一方程表明,向量 y 应由两个其他向量的线性组合来逼近。一个是元素全部为 1 的常向量,另一个是带有分量 exp(-t) 的向量。未知系数 c 1 c 2 可以通过执行最小二乘拟合来计算,这样会最大限度地减小数据与模型偏差的平方和。在两个未知系数的情况下有六个方程,用 6×2 矩阵表示。
E = [ones(size(t)) exp(-t)]
E = 6×2
1.0000 1.0000
1.0000 0.7408
1.0000 0.4493
1.0000 0.3329
1.0000 0.2019
1.0000 0.1003
        使用反斜杠运算符获取最小二乘解。
c = E\y
c = 2×1
0.4760
0.3413
        也就是说,对数据的最小二乘拟合为
 
        以下语句按固定间隔的 t 增量为模型求值,然后与原始数据一同绘制结果:
T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

        E*c y 不完全相等,但差值可能远小于原始数据中的测量误差。如果矩形矩阵 A 没有线性无关的列,则该矩阵秩亏。如果 A 秩亏,则 AX = B 的最小二乘解不唯一。如果A 秩亏,则 A\B 会发出警告,并生成一个最小二乘解。您可以使用 lsqminnorm 求在所有解中具有最小范数的解 X

4欠定方程组

        本例演示了欠定方程组的解不唯一的情况。欠定线性方程组包含的未知数比方程多。MATLAB 矩阵左除运算求基本最小二乘解,对于 m × n 系数矩阵,它最多有 m 个非零分量。
以下是一个简单的随机示例:
R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)
R =
6 8 7 3
3 5 4 1
b =
7
8
线性方程组 Rp = b 有两个方程,四个未知数。由于系数矩阵包含较小的整数,因此适合使用 format 命令以有理格式显示解。通过以下命令可获取特定解
format rat
p = R\b
p =
0
17/7
0
-29/7
        其中一个非零分量为 p(2) ,因为 R(:,2) 是具有最大范数的 R 的列。另一个非零分量为 p(4) ,因为 R(:,4) 在消除 R(:,2) 后起控制作用。欠定方程组的完全通解可以通过 p 加上任意零空间向量线性组合来表示,可以使用 null 函数(使用请求有理基的选项)计算该空间向量。
Z = null(R,'r')
Z =
-1/2 -7/6
-1/2 1/2
1 0
0 1
        可以确认 R*Z 为零,并且残差 R*x - b 远远小于任一向量 x(其中x = p + Z*q)  由于 Z 的列是零空间向量,因此 Z*q 是以下向量的线性组合:
        为了说明这一点,选择任意 q 并构造 x,计算残差的范数。
q = [-2; 1];
x = p + Z*q;
format short
norm(R*x - b)
ans =
2.6645e-15
        如果有无限多个解,则最小范数解具有特别意义。您可以使用 lsqminnorm 计算最小范数最小二乘解。该解具有 norm(p) 的最小可能值。
p = lsqminnorm(R,b)
p =
-207/137
365/137
79/137
-424/137

5多右端线性方程组的求解

        某些问题涉及求解具有相同系数矩阵 A 但具有不同右端 b 的线性方程组。如果可以同时使用不同的 b 值,则可以将 b 构造为多列矩阵,并使用单个反斜杠命令求解所有方程组: X = A\[b1 b2 b3 …]。但是,有时不同的 b 值并非全部同时可用,也就是说,您需要连续求解若干方程组。如果使用斜杠 (/) 或反斜杠 (\) 求解其中一个方程组,则该运算符会对系数矩阵 A 进行分解,并使用此矩阵分解来求解。然而,随后每次使用不同的 b 求解类似方程组时,运算符都会对 A 进行同样的分解,而这是一次冗余计算。此问题的求解是预先计算 A 的分解,然后重新使用因子对 b 的不同值求解。但是,实际上,以这种方式预先计算分解可能很困难,因为需要知道要计算的分解(LU、LDL、Cholesky 等)以及如何乘以因子才能对问题求解。例如,使用 LU 分解,您需要求解两个线性方程组才能求解原始方程组 Ax = b:
[L,U] = lu(A);
x = U \ (L \ b);
        对于具有若干连续右端的线性方程组,建议使用 decomposition 对象求解。借助这些对象,您可利用预先计算矩阵分解带来的性能优势,而不必了解如何使用矩阵因子。您可以将先前的 LU 分解替换为:
dA = decomposition(A,'lu');
x = dA\b;
        如果您不确定要使用哪种分解,decomposition(A) 会根据 A 的属性选择正确的类型,类似于反斜杠的功能。
        以下简单测试验证了此方法可能带来的性能优势。该测试分别使用反斜杠 (\) 和 decomposition 对同一稀疏线性方程组求解 100 次。
n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);
% Backslash solution
tic
for k = 1:100
x = A\b;
end
toc
Elapsed time is 9.006156 seconds.
% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
x = dA\b;
end
toc
Elapsed time is 0.374347 seconds.
        对于这个问题,decomposition 求解比单独使用反斜杠要快得多,而语法仍然很简单。

6迭代法

        如果系数矩阵 A 很大并且是稀疏矩阵,分解方法一般情况下将不会有效。迭代方法可生成一系列近似解。MATLAB 提供了多个迭代方法来处理大型的稀疏输入矩阵。

7多线程计算

        对于许多线性代数函数和按元素的数值函数,MATLAB 软件支持多线程计算。这些函数将自动在多个线程上执行。要使函数或表达式在多个 CPU 上更快地执行,必须满足许多条件:
1 函数执行的运算可轻松划分为并发执行的多个部分。这些部分必须能够在进程之间几乎不通信的情况下执行。它们应需要很少的序列运算。
2数据大小足以使并发执行的任何优势在重要性方面超过对数据分区和管理各个执行线程所需的时间。例如,仅当数组包含数千个或以上的元素时,大多数函数才会加速。
3运算未与内存绑定;处理时间不受内存访问时间控制。一般而言,复杂函数比简单函数速度更快。如果启用多线程, inv lscov linsolve mldivide 将会对大型双精度数组(约 10,000 个元素或更多)大幅增加速度。

 

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/818986.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

测试|自动化测试(了解)

测试|自动化测试&#xff08;了解&#xff09; 1.什么是自动化测试☆☆☆☆ 自动化测试相当于把人工测试手段进行转换&#xff0c;让代码执行。 2.自动化测试的分类☆☆☆☆ 注&#xff1a;这里只是常见的自动化测试&#xff0c;并不全部罗列。 1.单元自动化测试 其中Java…

分布式开源监控Zabbix实战

Zabbix作为一个分布式开源监控软件&#xff0c;在传统的监控领域有着先天的优势&#xff0c;具备灵活的数据采集、自定义的告警策略、丰富的图表展示以及高可用性和扩展性。本文简要介绍Zabbix的特性、整体架构和工作流程&#xff0c;以及安装部署的过程&#xff0c;并结合实战…

数据结构 | Radix Tree 树

什么是基数树&#xff1f; 基数树是一种多叉搜索树&#xff0c;数据位于叶子节点上&#xff0c;每一个节点有固定的2^n个子节点&#xff08;n为划分的基大小&#xff0c;当n为1时&#xff0c;为二叉树&#xff09;。 什么为划分的基&#xff1f; 以一个64位的长整型为例&#x…

oracle 19c打补丁遭遇OPATCHAUTO-72043OPATCHAUTO-68061

最近&#xff0c;在AIX上的新装oracle 19C数据库基础版本&#xff0c;使用opatchauto打PSU补丁集35037840时遇到了OPATCHAUTO-72043报错&#xff0c;无法正常应用GI补丁。 一、环境描述 操作系统&#xff1a;AIX 数据库版本&#xff1a;oracle rac 19.3.0新装基础版 应用PS…

代码随想录第四十八天|198、213、337.打家劫舍

198.打家劫舍 你是一个专业的小偷&#xff0c;计划偷窃沿街的房屋。每间房内都藏有一定的现金&#xff0c;影响你偷窃的唯一制约因素就是相邻的房屋装有相互连通的防盗系统&#xff0c;如果两间相邻的房屋在同一晚上被小偷闯入&#xff0c;系统会自动报警。 给定一个代表每个…

【JAVASE】重载与递归

⭐ 作者&#xff1a;小胡_不糊涂 &#x1f331; 作者主页&#xff1a;小胡_不糊涂的个人主页 &#x1f4c0; 收录专栏&#xff1a;浅谈java &#x1f496; 持续更文&#xff0c;关注博主少走弯路&#xff0c;谢谢大家支持 &#x1f496; 重载与递归 1. 方法重载1.1 为什么需要方…

开源代码分享(9)—面向100%清洁能源的发输电系统扩展规划(附matlab代码)

1.背景介绍 1.1摘要 本文提出了一种新颖的建模框架和基于分解的解决策略&#xff0c;将随机规划&#xff08;SP&#xff09;和鲁棒优化&#xff08;RO&#xff09;相结合&#xff0c;以应对协调中长期电力系统规划中的多重不确定性。从独立系统运营商&#xff08;ISO&#xff…

CAD曲面建模

首先还是要在3d绘图的环境中进行 在网络菜单栏的左侧&#xff0c;有曲面建模相关的功能 包括旋转形成曲面、平移形成曲面&#xff0c;按照两条线形成网格&#xff0c;按照四条封闭的线形成曲面等等 首先演示旋转曲面 需要被旋转的对象&#xff0c;以及旋转轴 首先选择要旋转…

Leetcode | DP | 338. 198. 139.

338. Counting Bits 重点在于这张图。 从i1开始&#xff0c;dp的array如果i是2的1次方之前的数&#xff0c;是1 dp[i - 2 ^ 0]; 如果i是2的2次方之前的数&#xff0c;是1 dp[i - 2 ^ 1]; 如果i是2的3次方之前的数&#xff0c;是1 dp[i - 2 ^ 2]; 198. House Robber 如果…

Hadoop学习指南:探索大数据时代的重要组成——Hadoop运行模式(下)

Hadoop运行模式(下&#xff09; 前言2.6 配置历史服务器1&#xff09;配置mapred-site.xml2&#xff09;分发配置3&#xff09;在hadoop102启动历史服务器4&#xff09;查看历史服务器是否启动5&#xff09;查看JobHistory 2.7 配置日志的聚集1&#xff09;配置yarn-site.xml2&…

【ChatGPT辅助学Rust | 基础系列 | Rust初相识】Rust简介与环境配置

教程目录 前言一&#xff0c;Rust简介1&#xff0c;Rust的历史2&#xff0c;Rust的特性3&#xff0c;为什么选择Rust4&#xff0c;Rust可以做什么 二&#xff0c; Rust环境配置1&#xff0c;windows11安装2&#xff0c;Linux安装 三&#xff0c;安装IDE 前言 Rust是一种系统编…

智慧水务配电能效系统的开发与功能介绍

随着城市化进程的步伐大大变快&#xff0c;城市建设与科学信息技术的融合程度也在不断提升&#xff0c;尤其是大数据信息技术的迅猛发展&#xff0c;为民生工程由信息化向智能化转型提供了条件。以城市的水务系统为例&#xff0c;依托大数据信息技术构建智慧水务系统是智慧城市…

从 DejaVu 改为 Noto,Ubuntu 23.10 发行版计划调整字体包

近日消息&#xff0c;代号为“Mantic Minotaur”的 Ubuntu 23.10 发行版计划调整字体包&#xff0c;从 DejaVu 修改为 Noto。 近日消息&#xff0c;代号为“Mantic Minotaur”的 Ubuntu 23.10 发行版计划调整字体包&#xff0c;从 DejaVu 修改为 Noto。 Ubuntu 开发团队表示为…

数字人第一剑,先斩“尹天仇”

第一波被数字人抢走饭碗的人类&#xff0c;不是带货主播&#xff0c;也不是虚拟偶像&#xff0c;而是好莱坞的群演们。 过去几个月&#xff0c;数千名好莱坞演员罢工&#xff0c;并出现在Netflix奈飞、华纳兄弟、Discovery、亚马逊、派拉蒙和NBC环球的办公大楼示威。 此次好莱坞…

哈工大计算机网络课程网络安全基本原理详解之:密钥分发中心与公钥认证中心

哈工大计算机网络课程网络安全基本原理详解之&#xff1a;密钥分发中心与公钥认证中心 在介绍密钥分发中心的概念前&#xff0c;先来回顾一下之前介绍的身份认证协议AP4.0&#xff1a;利用随机数R来避免“回放攻击”&#xff0c;并借助于对称加密算法来保证R的加密传输和解密&…

用python编写一个小程序,如何用python编写软件

大家好&#xff0c;给大家分享一下用python编写一个小程序&#xff0c;很多人还不知道这一点。下面详细解释一下。现在让我们来看看&#xff01; 1、python可以写手机应用程序吗&#xff1f; 我想有人曲解意思了&#xff0c;人家说用python开发渣蔽一个手机app&#xff0c;不是…

“ARTS挑战:探索技术,分享思考“

文章目录 前言一、学习的内容二、遇到的困难及解决办法三、学习打卡成果展示四、学习技巧的总结五、未来学习打卡计划后记 关于 ARTS 的释义 ● Algorithm: 每周至少做一个 LeetCode 的算法题 ● Review: 阅读并点评至少一篇英文技术文章 ● Tips: 学习至少一个技术技巧 ● Sha…

大数据面试实际场景类问题怎么准备,实在遭不住了!

300万字&#xff01;全网最全大数据学习面试社区等你来&#xff01; 前几天我发了一个面试总结的文章&#xff0c;实际场景类问题在大数据领域面试的占比越来越多。大家面试开始「务实」&#xff0c;在实际场景问题中&#xff0c;顺便增加对框架原理性内容的考察&#xff0c;这…

Docker中gitlab以及gitlab-runner的安装与使用

1、本文主要讲述如何使用Docker安装gitlab以及gitlab-runner&#xff0c;并且会讲述gitlab-runner如何使用 2、gitlab部分不需要修改过多的配置即可使用&#xff0c;本文未讲述https配置&#xff0c;如有需求&#xff0c;可自行百度 3、Docker如何安装可以自行百度 一、Docker安…

react中的高阶组件理解与使用

一、什么是高阶组件&#xff1f; 其实就是一个函数&#xff0c;参数是一个组件&#xff0c;经过这个函数的处理返回一个功能增加的组件。 二、代码中如何使用 1&#xff0c;高级组件headerHoc 2&#xff0c;在普通组件header中引入高阶组件并导出高阶组件&#xff0c;参数是普…