IT数码 购物 网址 头条 软件 日历 阅读 图书馆
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
图片批量下载器
↓批量下载图片,美女图库↓
图片自动播放器
↓图片自动播放器↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁
 
   -> Python知识库 -> 迭代法求解非线性方程组(含python代码) -> 正文阅读

[Python知识库]迭代法求解非线性方程组(含python代码)

1. 迭代法求解非线性方程组的原理

????????参考西安交大数值分析教材

?

?

?

?

?

?

?

?

?

?

2.?迭代法求解非线性方程组的计算过程

????????牛顿法求解非线性方程组的计算过程如下

?????????弦割法与牛顿法类似,弦割法将牛顿法中的偏导数用差商替代,即将上述牛顿法中的第(ii)步和第(iii)步转化为以下迭代格式

?

?????????布罗伊登法求解非线性方程组的计算过程如下

?3. 程序使用说明

????????程序使用方法为:(1)定义自变量;(2)使用(1)中的自变量定义非线性方程组列表;(3)给定初值,精度以及迭代次数;(4)调用方程netwon(),broyden(),string_cut(),分别对应牛顿法,布罗伊登法和弦割法。具体如下所示

4. python实现代码

?

'''
使用牛顿法、弦割法、布罗伊登法求解非线性方程组
'''
import sympy
import pprint
import numpy as np

#计算向量范数
def Norm2(p):
    sum_of_p = sum([i ** 2 for i in p])
    return sum_of_p ** 0.5

#牛顿法
def newton(A, x_0, args, prec=0.00000001, n=10):
    funcs=sympy.Matrix(A)
    args=args
    jacobianMatrix=funcs.jacobian(args) #计算雅可比矩阵
    #pprint.pprint(jacobianMatrix.inv())

    x_0=sympy.Matrix(x_0)
    x = args

    x_pre = x_0
    for k in range(0, n): #迭代求解
        b = - funcs.subs(zip(x, x_pre)) #计算f(x【k】)
        J = jacobianMatrix.subs(zip(x, x_pre))  #计算雅可比矩阵Jf(x【k】)
        deltx = sympy.Matrix(J.inv() * b) #使用雅可比矩阵的逆求解deltx
        x_new = x_pre + deltx  #迭代跟新x
        if Norm2(funcs.subs(zip(x, x_new))) < prec: #如果满足精度要求则返回
            #pprint.pprint(x_new)
            #print(k)
            return([x_new,k])
        x_pre = x_new
    return("迭代次数较少,无法满足精度要求")

#弦割法
def string_cut(A, x_0, args, prec=0.00000001, n=10):
    funcs = sympy.Matrix(A)
    #args = args
    h = 0.1 #设置h的大小
    e = np.eye(len(A))*h
    #pprint.pprint(e)
    x_0 = sympy.Matrix(x_0)
    x = args

    x_pre = x_0
    c=[]
    for k in range(0,n):
        for p in range(len(A)):
            c.append([A[p].subs(zip(x, (x_pre + sympy.Matrix(e[j]))))
                      for j in range(0, len(A))])
        fij = sympy.Matrix(c) #计算fij【k】
        #pprint.pprint(fij)
        b = funcs.subs(zip(x, x_pre)) #计算fi【k】
        z = fij.inv() * b  #求解z
        deltx = h * z / (sum(z) - 1)  #计算deltx
        x_new = x_pre + deltx #更新x
        if Norm2(funcs.subs(zip(x, x_new))) < prec: #满足精度要求则返回
            #pprint.pprint(x_new)
            #print(k)
            return([x_new,k])
        x_pre = x_new
        c = []
    return("迭代次数较少,无法满足精度要求")

def broyden(A, x_0, args, prec=0.00000001, n=10):
    funcs = sympy.Matrix(A)
    args = args
    jacobianMatrix = funcs.jacobian(args)

    x_0 = sympy.Matrix(x_0)
    x = args

    b = funcs.subs(zip(x, x_0))
    J0 = jacobianMatrix.subs(zip(x, x_0))
    x_1 = x_0 - J0.inv() * b
    prec = prec
    #以下得到两个初始迭代值和初始A【0】矩阵
    x_pre = x_0
    x_now = x_1
    J = jacobianMatrix.subs(zip(x, x_pre)).inv() #A【0】的逆
    #迭代计算
    for k in range(0, n):
        s = x_now - x_pre #计算s【k】
        y = funcs.subs(zip(x, x_now)) - funcs.subs(zip(x, x_pre)) #计算y【k】

        temp1 = (s - J * y) * s.T * J
        temp2 = (s.T * J * y)
        num = temp2 ** (-1)
        e=np.eye(len(A))
        n_num=np.array(num)
        e_num=e*(n_num)
        diag = sympy.diag(e_num)
        temp = temp1 * diag  #将除以一个数变为乘其逆的对角单位矩阵
        J_new = J + temp #得到A【k】的逆
        x_new = x_now - J_new * funcs.subs(zip(x, x_now)) #更新x的值
        if Norm2(funcs.subs(zip(x, x_new))) < prec: #如果满足精度要求则返回
            #pprint.pprint(x_new)
            #print(k)
            return [x_new, k]
        x_pre = x_now
        x_now = x_new
        J = J_new
    return ("迭代次数较少,无法满足精度要求")


if __name__ == "__main__":
    #---------------------------计算实习7.3.1-------------------------------
    x1, x2, x3 = sympy.symbols("x1 x2 x3") #定义自变量
    args = sympy.Matrix([x1, x2, x3]) #定义自变量参数矩阵
    #定义非线性方程组列表
    A = [x1 ** 2 + x2 ** 2 + x3 ** 2 - 1.0,
         2 * x1 ** 2 + x2 ** 2 - 4 * x3 ** 2,
         3 * x1 ** 2 - 4 * x2 ** 2 + x3 ** 2]
    x_0 = [1.0, 1.0, 1.0] #设置初始向量
    prec = 0.00000001 #设置精度
    n=20 #设置迭代次数
    p = newton(A, x_0, args, prec, n) #牛顿法求解
    pprint.pprint(p)
    p = broyden(A, x_0, args, prec, n) #布罗伊登法求解
    pprint.pprint(p)
    p = string_cut(A, x_0, args, prec, n) #弦割法求解
    pprint.pprint(p)
    #----------------------------计算实习7.3.2-------------------------------
    # x1, x2 = sympy.symbols("x1 x2")
    # args = sympy.Matrix([x1, x2])
    # A=[sympy.cos(x1**2+0.4*x2)+x1**2+x2**2-1.6, 1.5*x1**2-1/0.36*x2**2-1.0]
    # x_0=[1.04, 0.47]
    # prec = 0.00000001
    # n=20
    # p = newton(A, x_0, args, prec, n)
    # pprint.pprint(p)
    # p = broyden(A, x_0, args, prec, n)
    # pprint.pprint(p)
    # p = string_cut(A, x_0, args, prec, n)
    # pprint.pprint(p)

?

  Python知识库 最新文章
Python中String模块
【Python】 14-CVS文件操作
python的panda库读写文件
使用Nordic的nrf52840实现蓝牙DFU过程
【Python学习记录】numpy数组用法整理
Python学习笔记
python字符串和列表
python如何从txt文件中解析出有效的数据
Python编程从入门到实践自学/3.1-3.2
python变量
上一篇文章      下一篇文章      查看所有文章
加:2022-03-21 20:45:38  更:2022-03-21 20:50:03 
 
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁

360图书馆 购物 三丰科技 阅读网 日历 万年历 2024年12日历 -2024/12/29 12:46:24-

图片自动播放器
↓图片自动播放器↓
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
图片批量下载器
↓批量下载图片,美女图库↓
  网站联系: qq:121756557 email:121756557@qq.com  IT数码
数据统计