Assignment1 LU分解
- codes
- 参考文献
Assignment2 迭代法
- Problem
- 数据生成
- norm
- slow_Jacobi
- 清真_Jacobi
- Gauss-Seidel
- SOR超松弛法
- 参考文献
Assignment3CGQR
- Problem
- Conjugate Gradient Method 共轭梯度法
- QR Method
- 参考文献

Assignment1: LU分解

老师让只交.py，

于是很多东西直接在注释里写了

codes

# -*-encoding:utf-8-*-
# Numerical Computation & Optimization
# homework1: LU decomposition
# cww97
# 2017/09/27
# from cww97jh@gmail.com
# to num_com_opt@163.com
from numpy import *
n = 100def lu(A):  #"""Doolittle decomposition(course book page45)I just copy the formula from the course book:param A: the Coefficient matrix:return: L, U the lower & upper matrix"""L = mat(eye(n))  # low matrixU = mat(zeros((n, n)))   # upper matrixfor k in range(n):  # k-th stepfor i in range(k, n):  # calc k-th row of UU[k, i] = A[k, i] - sum(L[k, j]*U[j, i] for j in range(k))for i in range(k+1, n):  # k-th col of LL[i, k] = (A[i, k] - sum(L[i, j]*U[j, k] for j in range(k))) / U[k, k]return L, Udef calc(A, B):"""calc the equation Ax = b(course book page46)I just copy the formula from the course book:param L: lower triangle matrix:param U: upper triangle matrix:param B: B = A * X:return: the x of A * X = B"""L, U = lu(A)#  first, calc L * Y = B, get YY = mat(zeros((n, 1)))for k in range(n):Y[k, 0] = B[k, 0] - sum(L[k, j]*Y[j, 0] for j in range(k))#  then, calc U * X = Y, get XX = mat(zeros((n, 1)))for k in range(n-1, -1, -1):X[k, 0] = (Y[k, 0] - sum(U[k, j]* X[j, 0] for j in range(k+1, n))) / U[k, k]return X"""
sample data on course book, for debugA = mat([[2, 1, 5],[4, 1, 12],[-2, -4, 5]])B = mat([[11], [27], [12]])
"""
if __name__ == '__main__':# Generate a random matrix M &M = mat(random.randint(1, 100, size=(n, n)))# A = M + I(Identity matrix)A = M + mat(eye(n))print('---------matrix A:----------\n', A,)# Generate a vector x = (1,2,··· ,100).TX = mat([[i+1] for i in range(n)])print('---------vector X:----------\n', X.T, '.T')# Generate the vector b as b = A * xB = A * Xprint('-------B = A * X:-----------\n', B.T, '.T')x = calc(A, B)print('calc the equation A * x = B, x:\n', x.T, '.T')"""
finished this, when I need to calc some equation,
I think I may more likely to use this:x = np.linalg.solve(A, b)However, If I use this, this homework may cannot pass>_<!
"""

参考文献

python 矩阵操作

Assignment2: 迭代法

TAGS: 数值计算

数值计算与优化指导老师：王祥丰

陈伟文 10152510217

Problem

数据生成

def get_data():# Generate a random matrix M &A = mat(zeros((N, N)))for i in range(N):A[i, i] = 2for i in range(N-1):A[i, i + 1] = -1A[i + 1, i] = -1print('---------matrix A:----------\n', A,)# Generate a vector b = (1,1,··· ,1).Tb = mat(ones((N, 1)))print('---------vector X:----------\n', b.T, '.T')return A, b

norm

vector norm

matrix norm

slow_Jacobi

Wiki

核心算法

根据最后一个公式，写出如下代码

def norm(x):ans = 0for t in x:ans += square(t[0, 0])return sqrt(ans)def jacobi(A, b):n = len(b)x0 = mat(zeros((n, 1)))x1 = mat(zeros((n, 1)))d = 1ti = 0while d > eps:for i in range(n):tmp = 0for j in range(n):if j != i:tmp += A[i, j] * x0[j, 0]x1[i, 0] = 1./A[i, i] * (b[i, 0] - tmp)#print(x1.T)x0 = x1d = 1.*norm(A * x0 - b) / norm(b)ti += 1print('time = %d, d = %lf' % (ti, d))return x1

跑了好久好久好久，，，我问主席跑了多久

为啥他能跑这么快

于是乎我瞄了一眼他的代码，发现，他没像我这样一个值一个值的计算，而是使用numpy自带的矩阵运算
好的，我傻了，又自造车轮，太好了

于是我决定把上面的函数命名为slow_jacobi，来纪念我这个zz的操作

不过好歹等了十几分钟还是有结果出来的

运行速度慢了，不过迭代次数少了，不行，这么慢的代码我不能忍，重写

清真_Jacobi

用这个式子

Xk+1=D−1(b−Rxk)X^{k+1} = D^{-1} (b - Rx^k)

def jacobi(A, b):  # 这次吸取教训，多用矩阵运算# X^(k+1) = D^−1 (b − R * x^k)n = len(b)x0 = mat(zeros((n, 1)))x1 = mat(zeros((n, 1)))D = mat(diag(diag(A).tolist()))R = A - Dnorm_b = linalg.norm(b)d = 1; ti = 0while d > eps:x1 = D.I * (b - R * x0)x0 = x1; ti += 1d = linalg.norm(A * x0 - b) / norm_bprint('time = %d, delta = %.8lf' % (ti, d))return x0

迭代次数多了，不过10秒内出解

不要自造车轮，不要自造车轮，不要自造车轮

Gauss-Seidel

又是一波紧张刺激的抄公式~~（偷懒直接贴ppt了）~~

核心迭代还是抄一下吧

xk+1=BGxk+fGx^{k+1} = B_G x^k + f_G

BG=(D−L)−1，fG=(D−L)−1bB_G = (D - L) ^ {-1}，f_G = (D-L)^{-1}b

这个迭代可以用x迭代x

def gauss_seidel(A, b):# $x ^ {k + 1} = B_G x ^ k + f_G$# $B_G = (D - L) ^ {-1}, f_G = (D - L) ^ {-1}b$n = len(b)D = mat(diag(diag(A).tolist()))U = mat(zeros((n, n))) - triu(A, 1)L = mat(zeros((n, n))) - tril(A, -1)bg = (D - L).I * Ufg = (D - L).I * bnorm_b = linalg.norm(b)x = mat(zeros((n, 1)))d = 1; ti = 0while d > eps:x = bg * x + fgti += 1d = linalg.norm(A * x - b) / norm_bprint('time = %d, delta = %.8lf' % (ti, d))return x

迭代次数少了很多

SOR超松弛法

抄ppt

xk+1=Bwxk+fwx^{k+1} = B_w x^k + f_w

Bw=(D−wL)−1[(1−w)D+wU]B_w = (D - wL)^{-1}[(1-w)D + wU]

fw=w(D−wL)−1bf_w = w(D - wL)^{-1}b

关于ω\omega的取值（ppt上写成了ww），这篇论文，下载之后发现一共只有四页

一个跟实验报告一样的论文，还是c语言实现的

互动百科这么说

看来这个ω\omega取值，在11之间取吧22，就是看脸

def sor(A, b, w):# $x^{k+1} = B_w x^k + f_w$# $B_w = (D - wL)^{-1}[(1-w)D + wU]$# $f_w = w(D - wL)^{-1}b$n = len(b)D = mat(diag(diag(A).tolist()))U = mat(zeros((n, n))) - triu(A, 1)L = mat(zeros((n, n))) - tril(A, -1)bw = (D - w * L).I * ((1-w)*D + w*U)fw = w * (D - w*L).I * bnorm_b = linalg.norm(b)x = mat(zeros((n, 1)))d = 1; ti = 0while d > eps:x = bw * x + fwti += 1d = linalg.norm(A * x - b) / norm_bprint('w = %.2lf, time = %d' % (w, ti))return x

这个时候我已经不关心是否能出正解了，肯定是正解

为了测试下ω\omega的取值对迭代次数的影响

我又写了如下循环

    w = 1.while w <= 2:w += 0.1x = sor(A, b, w)

output:

w = 1.10, time = 11597
w = 1.20, time = 9448
w = 1.30, time = 7630
w = 1.40, time = 6071
w = 1.50, time = 4719
w = 1.60, time = 3536
w = 1.70, time = 2489
w = 1.80, time = 1554
w = 1.90, time = 696

是越大越快吗，我又改循环：

    w = 1.8while w < 2:w += 0.01x = sor(A, b, w)

得到如下输出

w = 1.81, time = 1466
w = 1.82, time = 1378
w = 1.83, time = 1292
w = 1.84, time = 1205
w = 1.85, time = 1120
w = 1.86, time = 1035
w = 1.87, time = 950
w = 1.88, time = 866
w = 1.89, time = 781
w = 1.90, time = 696
w = 1.91, time = 609
w = 1.92, time = 520
w = 1.93, time = 423
w = 1.94, time = 293
w = 1.95, time = 303
w = 1.96, time = 404
w = 1.97, time = 505
w = 1.98, time = 808
w = 1.99, time = 1515

我觉得没有必要继续枚举了

我只能说对于本题，ω\omega取1.941.94的时候速度最快

参考文献

[1] Lecture-3课程ppt

[2] wiki_jacobi

[3] 数值计算——线性方程组的迭代法

[4] 胡枫 ,于福溪《超松弛迭代法中松弛因子 ω的选取方法》青海师范大学学报 2006.01.13 42-46

[5] 互动百科_松弛法

[6] numpy文档_上三角矩阵

[7] numpy文档_norm

Assignment3：CG&QR

陈伟文 10152510217

2017/10/19

Problem

Calculate the equation Ax = b through Conjugate Gradient Method
and QR Method respectively

Conjugate Gradient Method 共轭梯度法

wiki，中文wiki只有寥寥几句，英文的比较详细

贴ppt

观察了一下每次循环的逻辑：

先x，用到了上一步的x和p，
r,用到上一步的r和p
p，用到刚刚算出的r和上一步的p

很显然可以发现，x，r，p只需要一个变量就可以完成循环

def conjugate_gradient(A, b):n = len(b)norm_b = linalg.norm(b)# generate x0, r0, p0x = mat(zeros((n, 1)))r = b - A * x; p = r  # 1# here we god = 1; ti = 0while d > eps:  #2ap = A * p  # 3a = (float)(1.*(r.T * r) / (p.T * ap))x += a * p  # 4r -= a * ap  # 5p = r + (float)(1.*(r.T * ap) / (p.T * ap)) * p  #6# for countingti += 1d = linalg.norm(A * x - b) / norm_bprint('time = %d, d = %f' % (ti, d))return x

我已经不关心方程解出来的正确性了（肯定是对的），看迭代次数吧

QR Method

先贴公式

先算Q

然后算R

（tips：这里的R特别容易算错）

最后：RX=QTbRX = Q^Tb

def QR(A, b):n = len(b)# calc QR = mat(zeros((n, n)))R[0, 0] = linalg.norm(A[:, 0])Q = A.copy()Q[:, 0] = A[:, 0] / linalg.norm(A[:, 0])for j in range(1, n):for i in range(j):Q[:, j] -= float(A[:, j].T * Q[:, i]) * Q[:, i]R[j, j] = linalg.norm(Q[:, j])Q[:, j] /= linalg.norm(Q[:, j])# calc Rfor i in range(n):for j in range(i+1, n):R[i, j] = float(A[:, j].T * Q[:, i])# calc x from Rx = Q^T * bb = Q.T * bx = mat(zeros((n, 1)))for i in range(n-1, -1, -1):x[i] = (b[i] - R[i, i+1:] * x[i+1:]) / R[i, i]return x

输出取了整：

calc the equation A * x = B, x:[[   50.    99.   147.   194.   240.   285.   329.   372.   414.   455.495.   534.   572.   609.   645.   680.   714.   747.   779.   810.840.   869.   897.   924.   950.   975.   999.  1022.  1044.  1065.1085.  1104.  1122.  1139.  1155.  1170.  1184.  1197.  1209.  1220.1230.  1239.  1247.  1254.  1260.  1265.  1269.  1272.  1274.  1275.1275.  1274.  1272.  1269.  1265.  1260.  1254.  1247.  1239.  1230.1220.  1209.  1197.  1184.  1170.  1155.  1139.  1122.  1104.  1085.1065.  1044.  1022.   999.   975.   950.   924.   897.   869.   840.810.   779.   747.   714.   680.   645.   609.   572.   534.   495.455.   414.   372.   329.   285.   240.   194.   147.    99.    50.]] .T

参考文献

[1] numpyt.dot 计算矩阵内积

[2] dot常见error

[3] numpy线性代数

[4] numpy行列操作

【数值计算】花式解线性方程组相关推荐

数值计算（一）之解线性方程组（高斯消去法，列选主元消去法，全选主元消去法，杜立特尔分解，克洛特分解，乔里斯基分解）
解线性方程组即解一个多元一次方程组,例如目录消去法分解法消去法原理没有学过高级的解法也没关系,凭借我们初高中的知识足以解决这个问题这是一个多元一次方程组,拥有n个未知量,也有n方程我们 ...
数值计算——雅可比迭代法解线性方程组(附代码)
1.雅克比迭代法的计算过程: (1).取初始向量: (1) (2 ...
数值计算实验2 解线性方程组实验（1）
文章目录实验目的: 实验内容: 需要word文件请访问 http://daxs.top 站内搜索实验名称或者实验内容访问文章并且下载附件即可. 实验目的: 进一步熟练掌握用Jacobi迭代法和Gau ...
二、解线性方程组的直接方法
https://zhuanlan.zhihu.com/p/30485749 设 n n n阶线性方程组: { a 11 x 1 + a 12 x 2 + . . . + a 1 n x n = b 1 ...
直接法 matlab,解线性方程组直接方法matlab用法.doc
解线性方程组直接方法matlab用法在这章中我们要学习线性方程组的直接法,特别是适合用数学软件在计算机上求解的方法. 2.1 方程组的逆矩阵解法及其MATLAB程序 2.1.3 线性方程组有解的判定 ...
【数理知识】《数值分析》李庆扬老师-第6章-解线性方程组的迭代法
第5章回到目录第7章第6章-解线性方程组的迭代法 6.1 迭代法的基本概念 6.2 雅可比迭代法与高斯-塞德尔迭代法 6.3 超松弛迭代法 6.4 共轭梯度法 6.1 迭代法的基本概念 6.2 ...
【数理知识】《数值分析》李庆扬老师-第5章-解线性方程组的直接方法
第4章回到目录第6章第5章-解线性方程组的直接方法 5.1 引言与预备知识 5.2 高斯消去法 5.3 矩阵三角分解法 5.4 向量和矩阵的范数 5.5 误差分析 5.1 引言与预备知识 5.2 ...
matlab算线性方程解,MATLAB计算方法3解线性方程组计算解法.pptx
第三章线性方程组数值解法解线性方程组 §3.1 直接法一. Gauss 消去法设有消元: 用Matlab实现顺序Gauss消去法在Matlab程序编辑器中输入:function x=nagauss ...
matlab lu分解求线性方程组_计算方法（二）直接三角分解法解线性方程组
封面是WH2里春希在编辑部的上司麻理前辈,有一说一,这条线的第一次H有点恶趣味,不是很喜欢. 一:概述矩阵分解我学过的挺多种,比如极分解,谱分解,满秩分解,正交三角分解还有这里的直接三角分解大部分我 ...

【数值计算】花式解线性方程组

Assignment1: LU分解

codes

参考文献

Assignment2: 迭代法

Problem

数据生成

norm

slow_Jacobi

清真_Jacobi

Gauss-Seidel

SOR超松弛法

参考文献

Assignment3：CG&QR

Problem

Conjugate Gradient Method 共轭梯度法

QR Method

参考文献

【数值计算】花式解线性方程组相关推荐

最新文章

热门文章