数值分析解线性方程组的编程实现(Hilbert)

列主元高斯消去（Gauss）

def Gauss(a,b,n):x=[0 for x in range(n)]for k in range(n-1):
#         找到列元素中最大的项maxIndex=kfor m in range(maxIndex+1,n):if abs(a[m][k])>abs(a[maxIndex][k]):maxIndex=mif a[maxIndex][k]==0:print("无解")returnif maxIndex!=k:for i in range(n):temp=a[maxIndex][i]a[maxIndex][i]=a[k][i]a[k][i]=temp#         消元计算for i in range(k+1,n):lik=a[i][k]/a[k][k]for j in range(k,n):a[i][j]=a[i][j]-lik*a[k][j]b[i]=b[i]-lik*b[k]if a[n-1][n-1]==0:print("无解")x[n-1]=b[n-1]/a[n-1][n-1]for i in range(n-1,0,-1):print(i-1)sum1=0for j in range(i,n):sum1+=a[i-1][j]*x[j]x[i-1]=(b[i-1]-sum1)/a[i-1][i-1]return xprint(x)
A=[[4,-1,2],[-1,-5,1],[2,1,6]]
b=[1,2,3]
Gauss(A,b,len(b))

SOR迭代法的实现：

def SOR(A,b,x,e,N,w):xr=[i for i in x]
#     迭代次数为N次for k in range(1,N+1):R=0#         算出每一行的增量for i in range(len(x)):sum1=0for j in range(len(x)):sum1+=A[i][j]*xr[j];R1=w*(b[i]-sum1)/A[i][i]if abs(R1)>R:R=R1xr[i]=xr[i]+R1if(R<=e):for i in range(len(xr)):print(xr[i],end="  ")returnprint("已达到最大迭代次数")for i in range(len(xr)):print(xr[i],end="  ")A=[[4,-1,2],[-1,-5,1],[2,1,6]]
b=[1,2,3]
x=[0,0,0]
e=0.01
N=10
w=1
SOR(A,b,x,e,N,w)

下面对Hilbert矩阵的求解使用gauss GS Jacobi SOR方法进行求解

import numpy as np
from math import isnan
import time
import sys
# 生成n阶的Hilbert矩阵
def generate(n):a = [[0.0 for i in range(n)] for j in range(n)]for i in range(n):for j in range(n):a[i][j]=1.0/(i+j+1)return a
# 计算a矩阵和b向量相乘
def mul(a,b):c=[0 for i in range(len(b))]for i in range(len(b)):sum1=0for j in range(len(b)):sum1+=a[i][j]*b[j]c[i]=sum1return c# noinspection PyUnresolvedReferences,PyTypeChecker
def Gauss(n):#     生成Hilbert矩阵a =generate(n)#     假设解的分量全部为1，算出bx = [1 for i in range(n)]b = mul(a,x)#     初始解向量x0 = [0 for i in range(n)]#   启动计时器t1=time.time()#     生成上三角矩阵for k in range(n - 1):for i in range(k + 1, n):lik = a[i][k] / a[k][k]for j in range(k, n):a[i][j] = a[i][j] - a[k][j] * likb[i] = b[i] - b[k] * lik#      根据上三角矩阵算出解x0[n - 1] = b[n - 1] / a[n - 1][n - 1]for i in range(n - 1, 0, -1):sum = 0for j in range(i, n):sum += a[i - 1][j] * x0[i]x0[i - 1] = (b[i - 1] - sum) / a[i - 1][i - 1]t2=time.time()for i in range(n):print(x0[i])print("误差为" ,end=" ")print(max([abs(np.array(x0)[i] - x[i]) for i in range(n)]))print("时间为" ,end=" ")print(t2-t1," 秒")# noinspection PyUnresolvedReferences
def Jacobi(n, e=0.00001):#     生成Hilbert矩阵a =generate(n)#     假设解的分量全部为1，算出bx = [1 for i in range(n)]b = mul(a,x)#     初始解向量x0 = [0 for i in range(n)]#   启动计时器t1=time.time()#     算法迭代开始。m=0while(True):m+=1x1 = [x0[i] for i in range(n)]for i in range(n):if isnan(x1[i]):print("数值无法计算")for k in range(n):print(x1[k])returnfor i in range(n):sum1 = 0for j in range(n):sum1 += a[i][j] * x0[j]x1[i] = x0[i] + (b[i] - sum1) / a[i][i]#             R为Xk和X(k-1)的误差R = max([abs(x0[i] - x1[i]) for i in range(n)])#         如果精度<e则输出结果if (R < e):t2=time.time()e = max([abs(x1[i] - 1) for i in range(n)])print("精度为e=", e, " 的解为:")for i in range(n):print(x1[i])print("时间为",end="  ")print(t2 - t1, " 秒")print("迭代次数为：", m)returnx0 = [x1[i] for i in range(n)]# noinspection PyUnresolvedReferences
def GS(n, e=0.00001, w=1.0):#     生成Hilbert矩阵a =generate(n)#     假设解的分量全部为1，算出bx = [1 for i in range(n)]b = mul(a,x)#     初始解向量x0 = [0 for i in range(n)]#       启动计时器t1=time.time()#    迭代开始m=0while(True):R=0m+=1for i in range(n):sum1 = 0for j in range(n):sum1 += a[i][j] * x0[j]R1 = w * (b[i] - sum1) / a[i][i]x0[i] = x0[i] + R1#          R为本次迭代解和实际解[1,1,1....]的误差if abs(R1)>R:R=R1#         如果精度<e则输出结果if (abs(R) < e):if(w==1.0):t2=time.time()e=max([abs(x0[i]-1) for i in range(n)])print("精度为e=", e, " 的解为:")for i in range(n):print(x0[i])print("时间为",end=" ")print(t2 - t1, " 秒")print("迭代次数为：", m)t2=time.time()return m,x0,t2-t1# 松弛因子
def SOR(n, e=0.00001):# 迭代次数min=sys.maxsize# 矩阵迭代效果最好时对应的松弛因子b=0.0# 对应的解为xx=[0 for i in range(n)]# 对应的时间t=0k=0.01# 循环调用GS方法传入松弛因子，最终找到迭代次数# 最少的对应的松弛因子for i in range(199):a,x0,t1=GS(n,e,k)if(a<min):b=kx=x0min=at=t1k+=0.01print("松弛因子为：", b)e=max([abs(x[i]-1) for i in range(n)])print("精度为：",e," 的解为")for i in range(n):print(x[i])print("时间为",t," 秒")print("迭代次数为：",min)def fun():print("矩阵为六阶的情况：")print("Gauss:")Gauss(6)print("Jacobi:")Jacobi(6,0.00001)print("GS:")GS(6,0.00001)print("SOR:")SOR(6,0.000001)# print("\n\n\n\n")# print("矩阵为25阶的情况")# print("Gauss:")# Gauss(25)# print("Jacobi:")# Jacobi(25, 0.0001)# print("GS:")# GS(25, 0.00001)# print("SOR:")# SOR(25, 0.00001)fun()

数值分析解线性方程组的编程实现(Hilbert)相关推荐

数值分析解线性方程组的直接法(一)
数值分析第三章引言高斯消去法高斯消去法的基本思想 n元线性方程组的高随消去法高斯消去法的计算步骤列主元高斯消去法列主元高斯消去法的基本思想算法步骤列主元高斯消元法的代码实例引言对于 ...
数值分析—行主元消去法解线性方程组—FORTRAN程序
数值分析-行主元消去法解线性方程组-FORTRAN程序 program main implicit none real8,dimension( :,: ),allocatable::A real8,d ...
克拉默法则C语言编程,FORTRAN编程：克拉默法则解线性方程组
FORTRAN编程:克拉默法则解线性方程组摘要:求解线性方程组的方法多种多样,例如:追赶法.高斯消去法.迭代法等等.我们在线性代数中学习过用克拉默法则来求解线性方程组,它旨在计算出几个矩阵的行列式即 ...
【数理知识】《数值分析》李庆扬老师-第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 ...
超松弛迭代法解线性方程组c语言,超松弛迭代法解线性方程组.doc
PAGE PAGE 2 姓名:___________________________ 设计题目:超松弛迭代法解线性方程组专业: 摘要本文是在matlab环境下熟悉的运用计算机编程语言并结合超松弛变 ...
选主元的高斯-约旦（Gauss-Jordan）消元法解线性方程组/求逆矩阵
选主元的高斯-约当(Gauss-Jordan)消元法在很多地方都会用到,例如求一个矩阵的逆矩阵.解线性方程组(插一句:LM算法求解的一个步骤),等等.它的速度不是最快的,但是它非常稳定(来自网上的定义 ...
计算方法（三）平方根法及其改进解线性方程组
一:概述本篇文章介绍解线性方程组的平方根法及改进平方根法,适用范围为系数矩阵为正定Hermite矩阵(下称H阵)的线性方程组.这个方法的理论依据我觉得是来自Schur引理的H阵结构定理,从这个角度我 ...
选主元的高斯-约当（Gauss-Jordan）消元法解线性方程组和求逆矩阵
选主元的高斯-约当(Gauss-Jordan)消元法在很多地方都会用到,例如求一个矩阵的逆矩阵.解线性方程组(插一句:LM算法求解的一个步骤),等等.它的速度不是最快的,但是它非常稳定(来自网上的定义 ...
高斯消去法解线性方程组的fortran程序实现
高斯消去法解方程组的fortran程序实现许多实际问题的解决,常常要化为求解线性代数方程组:例如,用最小二乘法处理测量结果和用差分法求解偏微分方程时,都会得到线性方程组:同时,很多物理学的问题最后也 ...

数值分析解线性方程组的编程实现(Hilbert)

数值分析解线性方程组的编程实现(Hilbert)相关推荐

最新文章

热门文章

数值分析 解线性方程组的编程实现(Hilbert)

数值分析 解线性方程组的编程实现(Hilbert)相关推荐

最新文章

热门文章

数值分析解线性方程组的编程实现(Hilbert)

数值分析解线性方程组的编程实现(Hilbert)相关推荐