数值分析第五章节 用Python实现解线性方程组的直接解法

news2024/9/28 7:16:13

参考书籍:数值分析 第五版 李庆杨 王能超 易大义编 第5章 解线性方程组的直接解法
文章声明:如有发现错误,欢迎批评指正

文章目录

  • 引言与预备知识
  • 高斯消去法
    • 列主元消去法
  • 矩阵三角分解法
    • 杜利特尔分解法
    • 平方根法
  • 向量和矩阵的范数
  • 误差分析

引言与预备知识

5.1.1:关于线性方程组的数值解法一般有两类:直接法和迭代法。5.1.2向量和矩阵:非奇异矩阵:设 A ∈ R n × n A\in R^{n\times n} ARn×n B ∈ R n × n B\in R^{n\times n} BRn×n。如果 A B = B A = I AB=BA=I AB=BA=I,则称 B B B A A A的逆矩阵,记为 A − 1 A^{-1} A1,且 ( A − 1 ) T = ( A T ) − 1 (A^{-1})^T=(A^T)^{-1} (A1)T=(AT)1。如果 A − 1 A^{-1} A1存在,则称 A A A为非奇异矩阵。如果 A , B ∈ R n × n A,B\in R^{n\times n} ABRn×n均为非奇异矩阵,则 ( A B ) − 1 = B − 1 A − 1 (AB)^{-1}=B^{-1}A^{-1} (AB)1=B1A1。其他跳过。5.1.3矩阵的特征值与谱半径: A A A的全体特征值称为 A A A的谱,记作 σ ( A ) \sigma(A) σ(A),即 σ ( A ) = { λ 1 , λ 2 , … , λ n } \sigma(A)=\{\lambda_1,\lambda_2,\dots,\lambda_n\} σ(A)={λ1,λ2,,λn}。记 ρ ( A ) = max ⁡ 1 ≤ i ≤ n ∣ λ i ∣ \rho(A)=\max\limits_{1\leq i\leq n}|\lambda_i| ρ(A)=1inmaxλi,称为矩阵 A A A的谱半径。其他跳过。5.1.4特殊矩阵:直接跳过。都是高等代数知识,所以我们直接跳过。

高斯消去法

5.2.1高斯消去法:十分简单大家都会,注意一下使用条件。定理:约化的主元素 a i i ( i ) ≠ 0 ( i = 1 , 2 , … , k ) a_{ii}^{(i)}\neq0(i=1,2,\dots,k) aii(i)=0(i=1,2,,k)的充要条件是矩阵 A A A的顺序主子式 D i ≠ 0 ( i = 1 , 2 , … , k ) D_i\neq0(i=1,2,\dots,k) Di=0(i=1,2,,k)(挺重要的,挺好证的)。矩阵的初等行变换其实就是左乘一个矩阵。 L ( 1 ) A ( 0 ) = A ( 1 ) , … , L ( n ) A ( n − 1 ) = A ( n ) 。 L = L ( n − 1 ) × ⋯ × L ( 1 ) , U = A ( n ) L^{(1)}A^{(0)}=A^{(1)},\dots,L^{(n)}A^{(n-1)}=A^{(n)}。L=L^{(n-1)}\times\dots\times L^{(1)},U=A^{(n)} L(1)A(0)=A(1),,L(n)A(n1)=A(n)L=L(n1)××L(1)U=A(n) 5.2.3列主元法:我们使用列主元法进行演示。5.2.2矩阵的三角分解:定理(矩阵的LU分解):设 A A A n n n阶矩阵,如果 A A A的顺序主子式 D i ≠ 0 ( i = 0 , 1 , … , n − 1 ) D_i\neq0(i=0,1,\dots,n-1) Di=0(i=0,1,,n1),则 A A A可分解为一个单位下三角矩阵 L L L和一个上三角矩阵 U U U的乘积,且这种分解是唯一的。5.2.1,5.2.3,5.2.2的顺序更有条理。

列主元消去法

我感觉我写得挺好,可以算作通用代码,前提必须保证有解。caozuo1是交换,caozuo2是归0。你把caozuo1删除了,就是高斯消去。这是第五章节例4,书上每步全都采用4位近似,我的更加准确一点。

def caozuo1(lt,x):
    global m
    a,b=x,abs(lt[x][x])
    for i in range(x+1,m):
        if abs(lt[i][x])>b:
            a,b=i,abs(lt[i][x])
    lt[x],lt[a]=lt[a],lt[x]
    return lt
def caozuo2(lt,x):
    global m,n
    for i in range(x+1,m):
        k=-lt[i][x]/lt[x][x];lt[i][x]=0
        for j in range(x+1,n):
            lt[i][j]+=lt[x][j]*k
    return lt
def solve(lt):
    global m,n
    x=[]
    for i in range(m-1,-1,-1):
        cnt=lt[i][n-1]
        for j in range(len(x)):
            cnt-=lt[i][n-j-2]*x[j]
        x.append(cnt/lt[i][i])
    return x[::-1]
m,n=map(eval,input().strip().split())
lt=[]
for _ in range(m):
    lt.append([eval(_) for _ in input().strip().split()])
for i in range(m-1):
    lt=caozuo1(lt,i)
    lt=caozuo2(lt,i)
print(solve(lt))

在这里插入图片描述

矩阵三角分解法

5.3.1直接三角分解法:1.不选主元的三角分解法。(杜利特尔分解法)2.选主元的三角分解法。我们使用杜利特尔分解法进行演示。
我感觉我写得挺好,可以算作通用代码,前提必须保证有解。这是第五章节例5,与书上的一摸一样。

杜利特尔分解法

def L(lt1,lt2,x):
    for i in range(x+1,n):
        cnt=0
        for j in range(x):
            cnt+=lt1[i][j]*lt2[j][x]
        lt1[i][x]=(lt[i][x]-cnt)/lt2[x][x]
    return lt1
def U(lt1,lt2,x):
    global lt,n
    for i in range(x,n):
        cnt=0
        for j in range(x):
            cnt+=lt1[x][j]*lt2[j][i]
        lt2[x][i]=lt[x][i]-cnt
    return lt2
n=int(input("系数方阵的行列数"))
lt1=[[0 for _ in range(n)] for _ in range(n)]
for _ in range(n):
    lt1[_][_]=1
lt2=[[0 for _ in range(n)] for _ in range(n)]
lt=[]
for _ in range(n):
    lt.append([eval(_) for _ in input().strip().split()])
for i in range(n-1):
    lt2=U(lt1,lt2,i)
    lt1=L(lt1,lt2,i)
lt2=U(lt1,lt2,n-1)
print("===============L===============")
for i in lt1:
    for j in i:
        print("{:>8.2f}".format(j),end="")
    print()
print("===============U===============")
for i in lt2:
    for j in i:
        print("{:>8.2f}".format(j),end=" ")
    print()

在这里插入图片描述
我们使用之前写的列主元消去法求解。y,x如下。
y:
在这里插入图片描述
x:
在这里插入图片描述
5.3.2平方根法。定理:设 A A A n n n阶对称矩阵,且 A A A的所有顺序主子式均不为零,则 A A A可分解为 A = L D L T A=LDL^T A=LDLT,其中 A A A为单位下三角矩阵, D D D为对角矩阵。定理:如果 A A A n n n阶对称正定矩阵,则存在一个实的非奇异下三角矩阵 L L L,使 A = L L T A=LL^T A=LLT,当限定 L L L的对角元素为正时,这种分解是唯一的。我们对它进行演示。

平方根法

自己随便捏的对称正定矩阵。代码变量LT是L的转置意思。

def caozuo(LT,x):
    global lt,n
    cnt=0
    for i in range(x):
        cnt+=pow(LT[i][x],2)
    LT[x][x]=pow(lt[x][x]-cnt,0.5)
    for i in range(x+1,n):
        cnt=0
        for j in range(x):
            cnt+=LT[j][x]*LT[j][i]
        LT[x][i]=(lt[x][i]-cnt)/LT[x][x]
    return LT
n=int(input("系数矩阵的行列数"))
lt=[]
for _ in range(n):
    lt.append([eval(_) for _ in input().strip().split()])
LT=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
    LT=caozuo(LT,i)
for i in LT:
    for j in i:
        print("{:>8.2f}".format(j),end=" ")
    print()

在这里插入图片描述
5.3.3 追赶法:追赶法更简单,我们就不演示。

向量和矩阵的范数

5.4.1向量范数:向量的 p p p范数: ∣ ∣ x ∣ ∣ p = ( ∑ i = 1 n ∣ x i ∣ p ) 1 p , p ∈ [ 1 , ∞ ) ||x||_p=(\sum\limits_{i=1}^n|x_i|^p)^{\frac{1}{p}},p\in[1,\infty) ∣∣xp=(i=1nxip)p1p[1,)。其他东西我们不管。5.4.2 矩阵范数: ∣ ∣ A ∣ ∣ ∞ = max ⁡ 1 ≤ i ≤ n ∑ j = 1 n ∣ a i j ∣ ||A||_{\infty}=\max\limits_{1\leq i\leq n}\sum\limits_{j=1}^n|a_{ij}| ∣∣A=1inmaxj=1naij(行范数)。 ∣ ∣ A ∣ ∣ 1 = max ⁡ 1 ≤ j ≤ n ∑ i = 1 n ∣ a i j ∣ ||A||_1=\max\limits_{1\leq j\leq n}\sum\limits_{i=1}^n|a_{ij}| ∣∣A1=1jnmaxi=1naij(列范数)。 ∣ ∣ A ∣ ∣ 2 = λ max ⁡ ( A T A ) ||A||_2=\sqrt{\lambda_{\max}(A^TA)} ∣∣A2=λmax(ATA) 。其他东西我们不管。求范数的程序简单所以我们不作演示。

误差分析

定义:设 A A A为非奇异矩阵,称数 c o n d ( A ) v = ∣ ∣ A − 1 ∣ ∣ v ∣ ∣ A ∣ ∣ v ( v = 1 , 2 , ∞ ) cond(A)_v=||A^{-1}||_v||A||_v(v=1,2,\infty) cond(A)v=∣∣A1v∣∣Av(v=1,2,)为矩阵A的条件数。 c o n d ( A ) 2 = ∣ ∣ A − 1 ∣ ∣ 2 ∣ ∣ A ∣ ∣ 2 = λ max ⁡ ( A T A ) λ min ⁡ ( A T A ) cond(A)_2=||A^{-1}||_2||A||_2=\sqrt{\frac{\lambda_{\max}(A^TA)}{\lambda_{\min}(A^TA)}} cond(A)2=∣∣A12∣∣A2=λmin(ATA)λmax(ATA) A A A为对称矩阵时 c o n d ( A ) 2 = ∣ λ 1 ∣ ∣ λ 2 ∣ cond(A)_2=\frac{|\lambda_1|}{|\lambda_2|} cond(A)2=λ2λ1 c o n d ( A ) ∞ = ∣ ∣ A − 1 ∣ ∣ ∞ ∣ ∣ A ∣ ∣ ∞ cond(A)_{\infty}=||A^{-1}||_{\infty}||A||_{\infty} cond(A)=∣∣A1∣∣A c o n d ( A ) cond(A) cond(A)越大方程越病态。

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

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

相关文章

剑指offer42.连续子数组的最大和

这道题挺简单的,看完题脑子里出现的想法就是用一个sum来把数组从前往后加,如果sum小于0,那么对于和来说是会减小的,所以这个时候直接把sum归零,然后从这个位置再往后加,用一个max_sum来记录sum的最大值&…

全面适配 | 走近openGauss数据库+鲲鹏欧拉操作系统

引入 全面适配 | openEuler操作系统 openGauss数据库 开篇 1、openEuler欧拉操作系统 百度百科:openEuler是覆盖全场景的创新平台,在引领内核创新,夯实云化基座的基础上,面向计算架构互联总线、存储介质发展新趋势,…

【LeetCode】剑指 Offer Ⅱ 第1章:整数(5道题) -- Java Version

题库链接:https://leetcode.cn/problem-list/e8X3pBZi/ 题目解决方案剑指 Offer II 001. 整数除法快速除 ⭐剑指 Offer II 002. 二进制加法模拟:StringBuilder ⭐剑指 Offer II 003. 前 n 个数字二进制中 1 的个数动规:res[i] res[i & (…

路由的配置

1、在router中设置路由导航跳转函数,在index.js文件中写这句话: 1.1 只要发生跳转, 就会调用这个函数: 1.2 导航的声明函数 2、访问系统访问控制系统如何形成 3、来一个导航守卫的案例:看看导航守卫的案例,写一个Main.Vue 和login…

【微软知识】微软相关技术知识分享

微软技术领域 一、微软操作系统: 微软的操作系统主要是 Windows 系列,包括 Windows 10、Windows Server 等。了解 Windows 操作系统的基本使用、配置和故障排除是非常重要的。微软操作系统(Microsoft System)是美国微软开发的Wi…

中文LLaMA模型和指令精调的Alpaca大模型:中文数据进行二次预训练,进一步提升了中文基础语义理解能力

项目设计集合(人工智能方向):助力新人快速实战掌握技能、自主完成项目设计升级,提升自身的硬实力(不仅限NLP、知识图谱、计算机视觉等领域):汇总有意义的项目设计集合,助力新人快速实…

Pytorch个人学习记录总结 10

目录 优化器 优化器 官方文档地址:torch.optimhttps://pytorch.org/docs/stable/optim.html Debug过程中查看的grad所在的位置: model --> Protected Atributes --> _modules --> ‘model’ --> Protected Atributes --> _modules -…

Linux学成之路(基础篇0(二十三)MySQL服务(主从MySQL服务和读写分离——补充)

目录 一、MySQL Replication概述 优点 异步复制(Asynchronous repication) 全同步复制(Fully synchronous replication) 半同步复制(Semisynchronous replication) 三、MySQL支持的复制 四、部署主从…

解决AttributeError: ‘DataParallel‘ object has no attribute ‘xxxx_fc1‘

问题描述 训练模型时,分阶段训练,第二阶段加载第一阶段训练好的模型的参数,接着训练 第一阶段训练,含有代码 if (train_on_gpu):if torch.cuda.device_count() > 1:net nn.DataParallel(net)net net.to(device)第二阶段训练…

WIZnet W5100S-EVB-Pico DHCP 配置教程(三)

DHCP协议介绍 什么是DHCP? 动态主机配置协议DHCP(Dynamic Host Configuration Protocol)是一种网络管理协议,用于集中对用户IP地址进行动态管理和配置。 DHCP于1993年10月成为标准协议,其前身是BOOTP协议。DHCP协议由…

【计算机网络】第 4 课 - 物理层

欢迎来到博主 Apeiron 的博客,祝您旅程愉快 ! 时止则止,时行则行。动静不失其时,其道光明。 目录 1、物理层的基本概念 2、物理层协议的主要任务 3、物理层任务 4、总结 1、物理层的基本概念 在计算机网络中,用来…

❤️创意网页:创意动态画布~缤纷移动涂鸦~图片彩色打码

✨博主:命运之光 🌸专栏:Python星辰秘典 🐳专栏:web开发(简单好用又好看) ❤️专栏:Java经典程序设计 ☀️博主的其他文章:点击进入博主的主页 前言:欢迎踏入…

【C/C++】#include<xxx.h>和#include“xxx.h“

2023年7月29日&#xff0c;周六晚上 今天下午和晚上花了不少时间去研究这个C/C的头文件以及#include<xxx.h>和#include"xxx.h"之间的区别&#xff0c;收获到了很多的很有用的知识。非常值得花时间来以博客的形式总结这些学习成果。 说实话&#xff0c;我挺想…

使用WGCLOUD监测安卓(Android)设备的运行状态

WGCLOUD是一款开源运维监控软件&#xff0c;除了能监控各种服务器、主机、进程应用、端口、接口、docker容器、日志、数据等资源 WGCLOUD还可以监测安卓设备&#xff0c;比如安卓手机、安卓设备等 我们只要下载对应的安卓客户端&#xff0c;部署运行即可&#xff0c;如下是下…

【Python】数据分析+数据挖掘——探索Pandas中的数据筛选

1. 前言 当涉及数据处理和分析时&#xff0c;Pandas是Python编程语言中最强大、灵活且广泛使用的工具之一。Pandas提供了丰富的功能和方法&#xff0c;使得数据的选择、筛选和处理变得简单而高效。在本博客中&#xff0c;我们将重点介绍Pandas中数据筛选的关键知识点&#xff…

x86架构ubuntu22用docker部署zsnes

0. 环境 x86 ubuntu22 1. 安装docker $ sudo apt remove docker docker-engine docker $ sudo apt update $ sudo apt install -y apt-transport-https ca-certificates curl software-properties-common$ curl -fsSL http://mirrors.aliyun.com/docker-ce/linux/ubuntu/gpg …

HiveSQL SparkSQL中常用知识点记录

目录 0. 相关文章链接 1. hive中多表full join主键重复问题 2. Hive中选出最新一个分区中新增和变化的数据 3. Hive中使用sort_array函数解决collet_list列表排序混乱问题 4. SQL中对小数位数很多的数值转换成文本的时候不使用科学计数法 5. HiveSQL & SparkSQL中炸裂…

leetcode 面试题 01.03. URL化

⭐️ 题目描述 &#x1f31f; leetcode链接&#xff1a;面试题 01.03. URL化 思路&#xff1a; 计算出空格的个数&#xff0c;我们可以知道最后一个字符的位置 endPos&#xff0c;再从后 end 向前遍历若不是空格正常拷贝&#xff0c;是空格则替换成 %20&#xff0c;最终当空格…

Linux系统编程之进程控制(上)

一、进程标识 1.pid 每个进程都有非负整数表示的唯一进程ID&#xff0c;即pid&#xff0c;其类型为pid_t类型。可用ps命令查看当前所有进程的信息&#xff0c;该命令可以加选项&#xff0c;一般使用ps -ef或ps axf(打印进程树)&#xff0c;查看当前系统所有进程的信息。需要注…

【Rust教程 | 基础系列 | Rust初相识】Rust简介与环境配置

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