乔列斯基(Cholesky)法解方程(python,数值积分)
第四课 乔列斯基(Cholesky)法解方程
首先要清楚二次型和正定矩阵
“二次型”可以定义为n个变量的二次表达式
如果这个二次型的所有变量X的值都等于或大于零，那么这个二次型就是“正的”。如果x1 = x2 = = xn =0的值为零时的正型称为“正定的”。非正定的正二次型称为“半正定”。
使用我们通常的向量和矩阵表达
二次型可以简便地表达为
其中[A]是“二次型Q(x)的矩阵”。如果|A|为零或非零，Q(x)是“奇异的”或非奇异的。
如果二次型{x}T [A]{x}是正定的，行列式
必须是正的
举例
三个值全是正的
因此二次型形式为
因为只有当x1 = x2 = x3 = 0时，它才能为零，所以它是正定的。
乔列斯基法
当系数矩阵[A]对称正定时，通过让[U]为[L]的转置，可以得到一个不同的因式分解。这个因式分解是可以写成
可以很容易得出，L11的平方=A11，L11L21=A12,L31L11=A13,所以，下面这个例子为
注意:如果对称矩阵不是正定的，那么乔列斯基法将不再适用，因为会分解出负值。
举例如下：
系数矩阵是正定对称的，则可以分解为
解出l值为
‘向前迭代法’
得到
‘向后迭代法’
得到
带宽方程
在许多工程应用中，当系数具有“带状”结构时。这意味着非零系数聚集在从矩阵左上角到右下角的对角线周围。一个典型的例子如图所示，在任何一行的“主对角线”的两边都不超过两个非零系数。这个系统的“带宽”是5。如果系数是对称的，则只需要存储和操作主对角线和每行两个以上的系数。在这种情况下，“半带宽”被称为2(不包括主对角线)
如果[A]是对称的，并且我们希望只在下面的三角形中存储非零项，如图2.1中的粗体所示，则只涉及15项，而不是36项(如果我们存储整个矩阵)。
这仍然有点低效，因为我们必须存储18个术语，包括前两行中不需要的零。然而，这种存储方法的优点是每行有三个项(半带宽加1)，这使得编程相当容易。主对角线位于第三列。
代码如下：
分为一个主程序和两个子程序，子程序为带宽形式的下三角分解cholin，和乔列斯基下的下三角向后迭代法chobac。
#使用带宽储存的乔列斯基LLT分解
import numpy as np
import math
import B
n=3
iw=2
iwp1=iw+1
lb=np.array([[0,0,16],[0,4,5],[8,-4,22]],dtype=np.float)
b=np.array([[4],[2],[5]],dtype=np.float)
print('带宽系数')
for i in range(1,n+1):print(lb[i-1,:])
print('右手边向量',b)
B.cholin(lb)
print('带宽形式下三角')
for i in range(1,n+1):print(lb[i-1,:])
B.chobac(lb,b)
print('解向量',b) cholin
def cholin(kb):n=kb.shape[0]iw=kb.shape[1]-1for i in range(1,n+1):x=0for j in range(1,iw+1):x=x+kb[i-1,j-1]**2kb[i-1,iw]=(kb[i-1,iw]-x)**0.5for k in range(1,iw+1):x=0if i+k<=n: if k!=iw: for l in range(iw-k,0,-1):x=x+kb[i+k-1,l-1]*kb[i-1,l+k-1]ia=i+kib=iw-k+1kb[ia-1,ib-1]=(kb[ia-1,ib-1]-x)/kb[i-1,iw]chobac
def chobac(kb,loads):n=kb.shape[0] iw=kb.shape[1]-1loads[0,0]=loads[0,0]/kb[0,iw]for i in range(2,n+1):x=0k=1if i<=iw+1:k=iw-i+2for j in range(k,iw+1):x=x+kb[i-1,j-1]*loads[i+j-iw-2,0]loads[i-1,0]=(loads[i-1,0]-x)/kb[i-1,iw]loads[n-1,0]=loads[n-1,0]/kb[n-1,iw]for i in range(n-1,0,-1):x=0l=i+iwif i>n-iw:l=nm=i+1for j in range(m,l+1):x=x+kb[j-1,iw+i-j]*loads[j-1,0]loads[i-1,0]=(loads[i-1,0]-x)/kb[i-1,iw]终端输出结果：
程序结果与计算结果一致
下面附带一个面向对象的cholesky分解
import numpy as np
import mathclass LinerSolver:def __init__(self, A, b):self.A = Aself.b = bdef CholeskiSolver(self):n = len(self.A)# 分解 [A] = [L] [L^T]for k in range(n):self.A[k,k] = math.sqrt(self.A[k,k] - np.dot(self.A[k,0:k], self.A[k,0:k]))for i in range(k+1,n):self.A[i,k] = (self.A[i,k] - np.dot(self.A[i,0:k], self.A[k,0:k])) / self.A[k,k]for k in range(1,n): self.A[0:k,k] = 0.0# 求解 [L]{y} = {b} for k in range(n):self.b[k] = (self.b[k] - np.dot(self.A[k,0:k], self.b[0:k])) / self.A[k,k]# 求解 [L^T]{x} = {y} for k in range(n-1,-1,-1):self.b[k] = (self.b[k] - np.dot(self.A[k+1:n,k], self.b[k+1:n])) / self.A[k,k]return self.bA = np.array([ [ 4, -1, 1],[-1, 4.25, 2.75],[1, 2.75, 3.5] ])b = np.array([4, 6, 7.25])cls = LinerSolver(A, b) #创建一个求解器的实例clsx = cls.CholeskiSolver() #调用Choleski法求解
print(x)乔列斯基(Cholesky)法解方程(python,数值积分)相关推荐
- 乔列斯基分解法求线性方程组的MATLAB程序实现
编写的 乔列斯基分解算法的MATLAB 程序如下: 功能:LL分解法求线性方程组AX=b的解调用格式:[X,L]= SymPosl (A,b) 其中, A:线性方程组的系数矩阵: b:线性方程组的常数 ...
- [矩阵的三角分解系列四] 乔累斯基(Cholesky)分解公式
乔累斯基分解公式 简介 LLT分解 证明 具体解法 稳定性 LDLT分解 证明 具体解法 例子 LLT分解 LDLT分解 引用 矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,LDU分解,杜利 ...
- 平方根法 乔累斯基分解Cholesky_解线性方程组的直接解法
平方根法 乔累斯基分解Cholesky_解线性方程组的直接解法 标签:计算方法实验 #include <stdio.h> #include <math.h>const int ...
- 乔利斯基三角分解_解线性方程组的直接法4.1-2.ppt
您所在位置:网站首页 > 海量文档  > 高等教育 > 微积分 解线性方程组的直接法4.1-2.ppt24页 本文档一 ...
- Cholesky分解、乔列斯基分解
一.简介 1.1 定理 Cholesky分解法 又叫 平方根法,是一种分解 正定Hermite矩阵 (即 A=AH\boldsymbol A = \boldsymbol A^\mathrm HA=AH ...
- 参数随机场，随机参数生成python代码，基于乔列斯基分解中点法分解
本人用python较多,但是matlib使用很少,没装在这个软件,故采用python语言重现编写,里面没有真实的单元坐标,只能假设三个单元进行测试,如存在错误请及时和我联系,平附上对应matlib代码 ...
- 乔利斯基三角分解_杜利特尔及乔利斯基三角分解
/**********改进乔利斯基三角分解**********/ void Improved_Cholesky() { int i,j,k; float t; L=(float **)malloc(s ...
- MATLAB之楚列斯基分解法（九）
楚列斯基分解法 楚列斯基分解是专门针对对称正定矩阵的分解.设A=aij∈Rn×nA=a_{ij}\in R^{n\times n}A=aij∈Rn×n是对称正定矩阵,A=RTRA = R^TRA=R ...
- 数值计算（一）之解线性方程组（高斯消去法，列选主元消去法，全选主元消去法，杜立特尔分解，克洛特分解，乔里斯基分解）
解线性方程组即解一个多元一次方程组,例如 目录 消去法 分解法 消去法 原理 没有学过高级的解法也没关系,凭借我们初高中的知识足以解决这个问题 这是一个多元一次方程组,拥有n个未知量,也有n方程 我们 ...
最新文章
- 实时获取ccd图像_薄膜瑕疵在线检测系统0.1mm检测精度_实时在线检测
- 微软软件推送服务器,向 UWP 应用添加推送通知 - Azure Mobile Apps | Microsoft Docs
- linux 瘦客户机系统,2X ThinClientOS基于Linux的瘦客户端系统 | MOS86
- OpenSUSE 13.1 和 OpenSUSE 12.3 用户如何安装 Cinnamon 2.2 桌面
- 安装phpredis扩展以及phpRedisAdmin工具
- SpringBoot之Dubbox
- Linux AIO的新归宿：io_uring（介绍，系统调用）
- C++ 常见错误（01） —— error LNK1104: 无法打开文件“avcodec.lib”
- POJ1679 The Unique MST —— 次小生成树
- Python Excel 批量生成二维码
- iphone机型分辨率
- 低信噪比环境下GPS信号识别捕获技术
- python 基础系列（四） — Python中的面向对象
- Magicbook 2018开启TPM2.0
- SpringBoot的优点及缺点
- 网络爬虫js逆向解决网站登录RSA加密问题，不使用selenium如何实现登录，session维持登录状态请求爬取
- .net之微信企业号开发(二) 企业号人员身份认证与开发
- jQuery学习笔记总结
- 嵌入式Linux系统开发入门宝典（第2版）
- ndarray 与 array 的区别 关系