数值分析第六章节 用Python实现解线性方程组的迭代法

news2024/9/28 11:19:10

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

文章目录

  • 迭代法的基本概念
  • 雅可比迭代法与高斯-塞格尔迭代法
    • 雅可比迭代法
    • 高斯-塞格尔迭代法

迭代法的基本概念

6.1.1引言:定义:(1)对于给定的线性方程组 x = B x + f x=Bx+f x=Bx+f,用公式 x ( k + 1 ) = B x ( k ) + f x^{(k+1)}=Bx^{(k)}+f x(k+1)=Bx(k)+f逐步带入求近似解的方法称为迭代法(或称为一阶定常迭代法,这里 B B B k k k无关)(2)如果 lim ⁡ k → ∞ x ( k ) \lim\limits_{k\rightarrow\infty}x^{(k)} klimx(k)存在(记为 x ∗ x^* x),称此迭代法收敛,显然 x ∗ x^{*} x就是此方程组的解,否则称此迭代法发散。6.1.2:向量序列与矩阵序列的极限:给定线性方程组 x = B x + f x=Bx+f x=Bx+f及一阶定常迭代法 x ( k + 1 ) = B x ( k ) + f x^{(k+1)}=Bx^{(k)}+f x(k+1)=Bx(k)+f式,对任意选取初始向量 x ( 0 ) x^{(0)} x(0),迭代法 x ( k + 1 ) = B x ( k ) + f x^{(k+1)}=Bx^{(k)}+f x(k+1)=Bx(k)+f式收敛的充要条件是矩阵 B B B的谱半径 ρ ( B ) < 1 \rho(B)<1 ρ(B)<1。其他跳过。

雅可比迭代法与高斯-塞格尔迭代法

雅可比迭代法

{ x ( 0 ) x ( k + 1 ) = B x ( k ) + f , k = 0 , 1 , … , x ( 0 ) 为初始向量, B = D − 1 ( L + U ) , f = D − 1 b \left\{\begin{matrix}x^{(0)}\\x^{(k+1)}=Bx^{(k)}+f,k=0,1,\dots,\end{matrix}\right.x^{(0)}为初始向量,B=D^{-1}(L+U),f=D^{-1}b {x(0)x(k+1)=Bx(k)+f,k=0,1,,x(0)为初始向量,B=D1(L+U),f=D1b
我感觉我写得挺好,可以算作通用代码,前提必须保证收敛。输入:输入系数矩阵行数,系数矩阵,初始向量,迭代次数。输出:解的向量。命名十分规范,懂了理论不难看懂。

def func1(B,x):
    #不通用的矩阵乘法
    global n
    lt=[]
    for i in range(n):
        cnt=0
        for j in range(n):
            cnt+=B[i][j]*x[j]
        lt.append(cnt)
    return lt
def func2(Bx,f):
    #不通用的矩阵加法
    global n
    lt=[]
    for i in range(n):
        lt.append(Bx[i]+f[i])
    return lt
n=int(input())
lt=[]
for _ in range(n):
    lt.append([eval(_) for _ in input().strip().split()])
D_inv=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
    D_inv[i][i]=1/lt[i][i]
L_sum_U=[[0 for _ in range(n)] for _ in range(n)]
for i in range(1,n):
    for j in range(i):
        L_sum_U[i][j]=-lt[i][j]
for i in range(n-1):
    for j in range(i+1,n):
        L_sum_U[i][j]=-lt[i][j]
B=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
    for j in range(n):
        B[i][j]=L_sum_U[i][j]*D_inv[i][i]
f=[0 for _ in range(n)]
for i in range(n):
    f[i]=D_inv[i][i]*lt[i][-1]
x=[eval(_) for _ in input().strip().split()]
num=int(input())
for _ in range(1,num+1):
    x=func2(func1(B,x),f)
print(x)

用的例1,一模一样。
在这里插入图片描述

高斯-塞格尔迭代法

{ x ( 0 ) x ( k + 1 ) = B x ( k ) + f , k = 0 , 1 , … , x ( 0 ) 为初始向量, B = ( D − L ) − 1 U , f = ( D − L ) − 1 b \left\{\begin{matrix}x^{(0)}\\x^{(k+1)}=Bx^{(k)}+f,k=0,1,\dots,\end{matrix}\right.x^{(0)}为初始向量,B=(D-L)^{-1}U,f=(D-L)^{-1}b {x(0)x(k+1)=Bx(k)+f,k=0,1,,x(0)为初始向量,B=(DL)1U,f=(DL)1b
我感觉我写得挺好,可以算作通用代码,前提必须保证收敛。输入:输入系数矩阵行数,系数矩阵,初始向量,迭代次数。输出:解的向量。命名十分规范,懂了理论不难看懂。

def func1(lt1,lt2):
    #矩阵乘法
    a,b=len(lt1),len(lt2[0])
    lt=[[0 for _ in range(b)] for _ in range(a)]
    for i in range(a):
        for j in range(b):
            for p in range(len(lt1[0])):
                lt[i][j]+=lt1[i][p]*lt2[p][j]
    return lt
def func2(lt1,lt2):
    #不通用的矩阵加法
    global n
    lt=[]
    for i in range(n):
        lt.append([lt1[i][0]+lt2[i][0]])
    return lt
n=int(input())
lt=[]
for _ in range(n):
    lt.append([eval(_) for _ in input().strip().split()])
D=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
    D[i][i]=lt[i][i]
L=[[0 for _ in range(n)] for _ in range(n)]
for i in range(1,n):
    for j in range(i):
        L[i][j]=-lt[i][j]
U=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n-1):
    for j in range(i+1,n):
        U[i][j]=-lt[i][j]
D_minus_L=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
    for j in range(n):
        D_minus_L[i][j]=D[i][j]-L[i][j]
#这里涉及一个求解下三角阵的逆矩阵
D_minus_L_inv=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
    for j in range(i):
        cnt=0
        for k in range(i):
            cnt-=D_minus_L[i][k]*D_minus_L_inv[k][j]
        D_minus_L_inv[i][j]=cnt/D_minus_L[i][i]
    D_minus_L_inv[i][i]=1/D_minus_L[i][i]
B=func1(D_minus_L_inv,U)
f=func1(D_minus_L_inv,[[lt[_][-1]] for _ in range(n)])
x=[[eval(_)] for _ in input().strip().split()]
num=int(input())
for _ in range(1,num+1):
    x=func2(func1(B,x),f)
print(x)

用的例1,一模一样。
在这里插入图片描述
就这样吧,剩下方法,自己研究。

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

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

相关文章

Python 进阶(三):正则表达式(re 模块)

❤️ 博客主页&#xff1a;水滴技术 &#x1f338; 订阅专栏&#xff1a;Python 入门核心技术 &#x1f680; 支持水滴&#xff1a;点赞&#x1f44d; 收藏⭐ 留言&#x1f4ac; 文章目录 1. 导入re模块2. re模块中的常用函数2.1 re.search()2.2 re.findall()2.3 re.sub()2.4…

020 - STM32学习笔记 - Fatfs文件系统(二) - 移植与测试

020 - STM32学习笔记 - Fatfs文件系统&#xff08;二&#xff09; - 移植与测试 上节学习了FatFs文件系统的相关知识&#xff0c;这节内容继续学习在STM32上如何移植FatFs文件系统&#xff0c;并且实现文件的创建、读、写与删除等功能。 一、FatFs文件系统移植 移植还是在之…

uniapp 全局数据(globalData)的设置,获取,更改

globalData&#xff0c;这是一种简单的全局变量机制。这套机制在uni-app里也可以使用&#xff0c;并且全端通用 因为uniapp基本上都是将页面&#xff0c;或者页面中相同的部分&#xff0c;进行组件化&#xff0c;所以会存在父&#xff0c;子&#xff0c;&#xff08;子&#xf…

AI革命:揭开微软无与伦比的AI技术面纱

来源&#xff1a;猛兽财经 作者&#xff1a;猛兽财经 2023年7月25日&#xff0c;全球科技行业的领导者之一微软(MSFT)公布了其2023财年第四季度的财报。 除了举世闻名的Windows操作系统&#xff0c;微软还通过笔记本电脑、个人电脑和服务器等产品改变了世界&#xff0c;该公司…

ARM裸机-2

1、搞清楚各种版本号 1.1、ARM的型号命名问题 ARM7和ARMv7不是一回事。 Cortex-A9比Cortex-A7更先出来。 型号很乱&#xff0c;初学者容易分不清哪个是哪个&#xff0c;比较迷茫。 1.2、ARM的几种版本号 ARM内核版本号&#xff08;ARM卖给别人的核心版本号&#xff09; …

IntelliJ IDEA 2023.2 主要更新了什么?(图文版)

&#x1f337;&#x1f341; 博主猫头虎 带您 Go to New World.✨&#x1f341; &#x1f984; 博客首页——猫头虎的博客&#x1f390; &#x1f433;《面试题大全专栏》 文章图文并茂&#x1f995;生动形象&#x1f996;简单易学&#xff01;欢迎大家来踩踩~&#x1f33a; &a…

无敌的服务注册中心Spring CloudAlibaba Nacos不进来看一看吗?

1、Nacos概述 Nacos支持发现、配置和管理几乎所有类型的服务。其Key features&#xff1a;Service Discovery And Service Health Check、Dynamic configuration management、Dynamic DNS service、Service governance and metadata management 官网 2、Nacos安装运行 安装 …

LearnOpenGL_Day1

文章目录 前期准备下载GLFW下载GLAD 引入库文件生成窗口重要概念——双缓冲 &#xff08;double buffer&#xff09;代码实现 学习总结 前期准备 下载GLFW GLFW DOWNLOAD 解压后使用CMake编译至新创建的bulid文件夹下&#xff1a; 下载GLAD 引入库文件 创建好工程&#x…

git的ssh方式对接码云

一、环境准备&#xff1a; 1、git下载&#xff0c;360管家或是百度。 2、vs2022&#xff0c;百度下载。 二、配置git&#xff1a; 1、打开准备存放文件的文件夹&#xff0c;右键&#xff0c;选择“Git Bash here”&#xff0c;弹出命令窗口&#xff0c; 输入&#xff1a;ss…

平板光波导中传播常数β的matlab计算(验证了是对的)

参照的是导波光学_王建(清华大学)的公式(3-1-2、3-1-3)&#xff0c;算的参数是这本书的图3-3的。 function []PropagationConstantsMain() clear;clc;close all lambda01.55;%真空或空气中的入射波长&#xff0c;单位um k02*pi/lambda0; m3;%导模阶数(需要人为指定) n11.62;%芯…

入门redis你一定需要知道的命令

1、各种数据类型的特点 字符串(string)&#xff1a;普通字符串&#xff0c;Redis中最简单的数据类型 哈希(hash)&#xff1a;也叫散列&#xff0c;类似于Java中的HashMap结构 列表(list)&#xff1a;按照插入顺序排序&#xff0c;可以有重复元素&#xff0c;类似于Java中的Li…

回归预测 | MATLAB实现WOA-ELM鲸鱼算法优化极限学习机多输入单输出回归预测

回归预测 | MATLAB实现WOA-ELM鲸鱼算法优化极限学习机多输入单输出回归预测 目录 回归预测 | MATLAB实现WOA-ELM鲸鱼算法优化极限学习机多输入单输出回归预测效果一览基本介绍程序设计参考资料 效果一览 基本介绍 Matlab实现WOA-ELM鲸鱼算法优化极限学习机多输入回归预测&#…

流数据湖平台Apache Paimon(二)集成 Flink 引擎

文章目录 第2章 集成 Flink 引擎2.1 环境准备2.1.1 安装 Flink2.1.2 上传 jar 包2.1.3 启动 Hadoop2.1.4 启动 sql-client 2.2 Catalog2.2.1 文件系统2.2.2 Hive Catalog2.2.3 sql 初始化文件 2.3 DDL2.3.1 建表2.3.2 修改表 2.4 DML2.4.1 插入数据2.4.2 覆盖数据2.4.3 更新数据…

idea插件开发-自定义语言4-Syntax Highlighter

SyntaxHighlighter用于指定应如何突出显示特定范围的文本&#xff0c;ColorSettingPage可以定义颜色。 一、Syntax Highter 1、文本属性键&#xfeff; TextAttributesKey用于指定应如何突出显示特定范围的文本。不同类型的数据比如关键字、数字、字符串等如果要突出显示都需…

ReID网络(一):MGN网络

Start MGN 1. 序言 现代基于感知的信息中&#xff0c;视觉信息占了80~85%。基于视觉信息的处理和分析被应用到诸如安防、电力、汽车等领域。 以安防市场为例&#xff0c;早在2017年&#xff0c;行业咨询公司IHS Market&#xff0c;我国在公共和私人领域安装有摄像头约1.76亿…

《TCP IP网络编程》第十三章

第 13 章 多种 I/O 函数 13.1 send & recv 函数 Linux 中的 send & recv&#xff1a; send 函数定义&#xff1a; #include <sys/socket.h> ssize_t send(int sockfd, const void *buf, size_t nbytes, int flags); /* 成功时返回发送的字节数&#xff0c;失败…

36.悬浮板

悬浮板 html部分 <div class"container"><div class"square"></div> </div>css部分 *{margin: 0;padding: 0; } body{background-color: #111;height: 100vh;overflow: hidden;display: flex;justify-content: center;align-it…

layui框架学习(33:流加载模块)

Layui中的流加载模块flow主要支持信息流加载和图片懒加载两部分内容&#xff0c;前者是指动态加载后续内容&#xff0c;示例的话可以参考csdn个人博客主页&#xff0c;鼠标移动到页面底部时自动加载更多内容&#xff0c;而后者是指页面显示图片时才会延迟加载图片信息。   fl…

苍穹外卖-day08

苍穹外卖-day08 本项目学自黑马程序员的《苍穹外卖》项目&#xff0c;是瑞吉外卖的Plus版本 功能更多&#xff0c;更加丰富。 结合资料&#xff0c;和自己对学习过程中的一些看法和问题解决情况上传课件笔记 视频&#xff1a;https://www.bilibili.com/video/BV1TP411v7v6/?sp…

第17节 R语言分析:生物统计数据集 R 编码分析和绘图

生物统计数据集 R 编码分析和绘图 生物统计学,用于对给定文件 data.csv 中的医疗数据应用 R 编码,该文件是患者人口统计数据集,包含有关来自各种祖先谱系的个体的标准信息。 数据集特征解释 脚本 output= file("Output.txt") # File name of output log sink(o…