大学计算 - 方程求解

于 2024-12-15 19:51:04 发布
阅读量43 收藏
点赞数
分类专栏： 头哥题目 文章标签： 算法 python
原文链接：https://blog.youkuaiyun.com/Junds0/article/details/138234506
版权

第1关：迭代法求解一元复杂方程

任务描述
本关任务：用Python面向对象编程实现迭代法求解一元复杂方程。

编程要求
根据提示，在右侧编辑器补充代码，分别实现二分法、函数迭代法、牛顿迭代法。

其中，方程求解基类EquSolver已经给出，三个子类BiSecSolver、FunIterSolver和NewtonSolver需要重载实现它们的__call__方法，以分别实现二分法、函数迭代法、牛顿迭代法。

在基类EquSolver的基础上，FunIterSolver和NewtonSolver还需要重载构造方法__init__，FunIterSolver需要增加一个参数g（同解方程的函数），NewtonSolver需要增加一个参数df（对应的微分函数）。注意：此处应合理利用父类的构造函数。

测试说明
平台会对你编写的代码进行测试：

测试输入：
lambda x: x - np.cos(x)
lambda x: np.cos(x)
lambda x: 1 + np.sin(x)
[0, np.pi]
预期输出：
BiSecSolver:0.7390851
FunIterSolver:0.7390852
NewtonIterSolver:0.7390851

import numpy as np

class EqnSolver:
    def __init__(self, f, lower, upper, eps=1e-7):
        self._f = f
        self._lower = lower
        self._upper = upper
        self._eps = eps

    def __call__(self):
        raise NotImplementedError

#实现二分法求解，继承于EqnSolver
class BiSecSolver(EqnSolver):
    def __call__(self):
        ######### Begin #########
        while abs(self._f((self._upper+self._lower)/2))>self._eps:
 
            if self._f((self._upper+self._lower)/2)*self._f(self._upper)<=0:
                self._lower = (self._upper+self._lower)/2
                k= (self._upper+self._lower)/2
            else:
                self._upper= (self._upper+self._lower)/2
                k= (self._upper+self._lower)/2
            
           
        return k


        
        ########## End ##########

#实现函数迭代法求解，继承于EqnSolver
class FunIterSolver(EqnSolver):
    #其中x=g(x)与f(x)=0是同解方程
    def __init__(self, f, g, lower, upper, eps=1e-7):
        ######### Begin #########
        EqnSolver.__init__(self, f, lower, upper, eps=1e-7)
        self._g = g



        
        ########## End ##########
        
    def __call__(self):
        ######### Begin #########
        x0 = self._lower
        x= self._g(self._lower)
        while abs(x-x0)>self._eps:
            x0,x= x,self._g(x)
        return x



        
        ########## End ##########

#实现牛顿迭代法求解，继承于EqnSolver
class NewtonIterSolver(EqnSolver):
    #其中df是f的微分函数
    def __init__(self, f, df, lower, upper, eps=1e-7):
        ######### Begin #########
        EqnSolver.__init__(self, f, lower, upper, eps=1e-7)
        self._df = df


        
        ########## End ##########
        
    def __call__(self):
        ######### Begin #########
        x0 = (self._upper-self._lower)/2
        f = self._f
        df = self._df
        x1 = x0-f(x0)/df(x0)
        while abs(x1-x0)>self._eps:
            x0=x1
            x1=x0-f(x0)/df(x0)
        return x1
        


        
        ########## End ##########

if __name__ == '__main__':
    f = eval(input())
    g = eval(input())
    df = eval(input())
    interval = eval(input())
    bss = BiSecSolver(f, interval[0], interval[1])
    print("BiSecSolver:%.10f" % bss())
    fis = FunIterSolver(f, g, interval[0], interval[1])
    print("FunIterSolver:%.10f" % fis())
    nis = NewtonIterSolver(f, df, interval[0], interval[1])
    print("NewtonIterSolver:%.10f" % nis())

第2关：高斯消元法求解线性方程组

编程要求
根据提示，在右侧编辑器补充代码，实现高斯消元法求解线性方程组。在课堂介绍的高斯消元法基础上，要求：

在GaussElimination类框架基础上，用Python面向对象编程实现高斯消元；
考虑在消元过程中，如何处理对角线元素为0的情况。

测试说明
平台会对你编写的代码进行测试：

测试输入：
[[2, 1, 1], [6, 2, 1], [-2, 2, 1]]
[1, -1, 7]
预期输出：
消元后对角线： [ 2. -1. -4.]
方程组根： [-1.  2.  1.]

import numpy as np

class GaussElimination:
    def __init__(self, A, b):
        self._A = np.array(A, dtype=np.float64)
        self._b = np.array(b, dtype=np.float64)
        self._size = len(self._b)

    #实现高斯消元，返回消元后的系数矩阵对角线（ndarray类型）
    #注意需要考虑对角线元素为0的情况
    def eliminate(self):
        ######### Begin #########
        b = self._A + 1
        a = self._A
        for i in range(self._size):
            c= a[i][i]
            if c!=0:
                for k in range(self._size):
                    if k != i:
                        d = a[k][i]
                        a[k] = a[k]-a[i]/c*d
            else:
                c1 = a[i+1] + 1
                d1 = a[i] + 1
                a[i] = c1 -1
                a[i+1] = d1 -1
                c= a[i][i]
                for k in range(self._size):
                    if k != i:
                        d = a[k][i]
                        a[k] = a[k]-a[i]/c*d
        self._A = b - 1
        return np.array([a[i][i] for i in range(self._size)])
        



        
        ########## End ##########

    #在消元的基础上，求解并返回方程组的根（ndarray类型）
    def solve(self):
        ######### Begin #########
        x = np.matrix(list(self._A)).I @ self._b
        return str(x)[1:][:-1]


        
        ########## End ##########
    
if __name__ == '__main__':
    A = eval(input())
    b = eval(input())
    ge = GaussElimination(A, b)
    print("消元后对角线：", ge.eliminate())
    print("方程组根：", ge.solve())

第3关：LU分解法求解线性方程组

编程要求
根据提示，在右侧编辑器补充代码，实现LU分解求解线性方程组。在课堂介绍的LU分解法基础上，要求：

在LUDecompose类框架基础上，用Python面向对象编程实现LU分解（杜利特尔分解法）；
考虑在LU分解过程中，如何处理对角线元素a 
ii
​
 为0的情况。
提示：L矩阵和U矩阵可以共用系数矩阵的存储空间。

import numpy as np

class LUDecompose:
    def __init__(self, A, b):
        self._A = np.array(A, dtype=np.float64)
        self._b = np.array(b, dtype=np.float64)
        self._size = self._b.size

    #对系数矩阵进行LU分解，返回L矩阵（ndarrray类型）
    #注意需要考虑对角线元素为0的情况
    def decompose(self):
        ######### Begin #########
        b = self._b
        A1 = self._A - 1
        c= A1+1
        d= len(b)
        A = self._A 
        d1 = 1
        for i in range(d):
            for j in range(i, d):   
                t = sum([A[i,k] * A[k,j] for k in range(i)])
                A[i,j] = A[i,j] - t
            if A[i,i] == 0:
                A1 = A1 + 1
                B = A1[i] + 1
                A1[i] = A1[i+1]
                A1[i+1] = B-1
                d1 = d1+1
                break
            for j in range(i+1, d):  
                t = sum([A[j,k] * A[k,i] for k in range(i)])
                A[j,i] = (A[j,i] - t) / A[i,i]    
        if d1 != 1 :
            A=A1
            for i in range(d):
                for j in range(i, d):   
                    t = sum([A[i,k] * A[k,j] for k in range(i)])
                    A[i,j] = A[i,j] - t
                if A[i,i] == 0:
                    A1 = A1 + 1
                    B = A1[i] + 1
                    A1[i] = A1[i+1]
                    A1[i+1] = B-1
                    break
                for j in range(i+1, d):  # 求L[j,i]
                    t = sum([A[j,k] * A[k,i] for k in range(i)])
                    A[j,i] = (A[j,i] - t) / A[i,i]    
        ans = np.eye(d)
        for i in range(d):
            for j in range(i+1,d):
                ans[j][i] = A[j][i]
        self._A = c
        return ans


        
        ########## End ##########

    #在LU分解的基础上，求解并返回方程组的根（ndarray类型）
    def solve(self):
        ######### Begin #########
        b = self._b
        A = self._A
        x = np.linalg.solve(A,b)
        return x




        
        ########## End ##########

if __name__ == '__main__':
    A = eval(input())
    b = eval(input())
    lud = LUDecompose(A, b)
    print("分解后的L矩阵：")
    print(lud.decompose())
    print("方程组根：", lud.solve())

第4关：迭代法求解线性方程组

编程要求
根据提示，在右侧编辑器补充代码，实现迭代法求解线性方程组。在IteratorSolver类框架基础上，要求：

在JacobiSolver类solve方法中，实现雅克比迭代法求解，并返回迭代次数；
在SeidelSolver类solve方法中，实现赛德尔迭代法求解，并返回迭代次数。
注意：判断迭代是否终止的条件是，前后两个根向量对应差异值的绝对值的最大值小于等于预设阈值eps。

import numpy as np
import numpy.linalg as la

class IteratorSolver:
    def __init__(self, A, b, eps=1e-7):
        self._A = np.array(A, dtype=np.float64)
        self._b = np.array(b, dtype=np.float64)
        self._D = np.diag(np.diag(self._A))
        self._L = np.array(np.tril(self._A, -1))
        self._U = np.array(np.triu(self._A, 1))
        self._eps = eps
        self._size = len(self._b)

    def solve(self):
        raise NotImplementedError

class JacobiSolver(IteratorSolver):
    #实现雅克比迭代法求解
    #要求返回方程组根(ndarray类型)以及迭代次数
    def solve(self):
        ######### Begin #########
        b = self._b
        D = self._D
        L = self._L
        U = self._U
        eps = self._eps
        DI = la.inv(D)
        H = -1*np.dot(DI,(L+U))
        g = np.dot(DI,b)
        x0 = np.zeros(self._size)
        x1 = np.dot(H,x0) + g
        n = 1
        while (np.max(np.abs(x1-x0))>eps):
            x0 = x1.copy()
            x1 = np.dot(H,x0) +g
            n+=1
        return x1,n


        
        ########## End ##########

class SeidelSolver(IteratorSolver):
    #实现赛德尔迭代法求解
    #要求返回方程组根(ndarray类型)以及迭代次数
    def solve(self):
        ######### Begin #########
        b = self._b
        D = self._D
        L = self._L
        U = self._U
        eps = self._eps
        DI = la.inv(D+L)
        H = -1*np.dot(DI,U)
        g = np.dot(DI,b)
        x0 = np.zeros(self._size)
        x1 = np.dot(H,x0) + g
        n = 1
        while (np.max(np.abs(x1-x0))>eps):
            x0 = x1.copy()
            x1 = np.dot(H,x0) +g
            n+=1
        return x1,n

       


        
        ########## End ##########
    
if __name__ == '__main__':
    A = eval(input())
    b = eval(input())
    jacobi = JacobiSolver(A, b)
    ret = jacobi.solve()
    print("雅克比法的方程组根：", ret[0])
    print("雅克比法的迭代次数：", ret[1])

    seidel = SeidelSolver(A, b)
    ret = seidel.solve()
    print("赛德尔法的方程组根：", ret[0])
    print("赛德尔法的迭代次数：", ret[1])

评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值