1. Toy模板网首页
  2. > Toy博客

数值分析上机题（上）

10月前 作者：Lemon_Yam 分类：Toy博客 阅读(41) 违法举报

这篇具有很好参考价值的文章主要介绍了数值分析上机题（上）。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

第一章

一、问题

设 S N = ∑ j = 2 N 1 j 2 − 1 S_N = \sum_{j=2}^N \frac{1}{j^2 - 1} SN​=∑j=2N​j2−11​，其精确值为 1 2 ( 3 2 − 1 N − 1 N + 1 ) \frac{1}{2}(\frac{3}{2} - \frac{1}{N} - \frac{1}{N+1}) 21​(23​−N1​−N+11​)

  • 编制按从大到小的顺序 S N = 1 2 2 − 1 + 1 3 2 − 1 + ⋯ + 1 N 2 − 1 S_N = \frac{1}{2^2 - 1} + \frac{1}{3^2 - 1} + \cdots + \frac{1}{N^2 - 1} SN​=22−11​+32−11​+⋯+N2−11​，计算 S N S_N SN​ 的通用程序
  • 编制按从小到大的顺序 S N = 1 N 2 − 1 + 1 ( N − 1 ) 2 − 1 + ⋯ + 1 N 2 − 1 S_N = \frac{1}{N^2 - 1} + \frac{1}{(N-1)^2 - 1} + \cdots + \frac{1}{N^2 - 1} SN​=N2−11​+(N−1)2−11​+⋯+N2−11​，计算 S N S_N SN​ 的通用程序
  • 按两种顺序分别计算 S 1 0 2 S_{10^2} S102​， S 1 0 4 S_{10^4} S104​， S 1 0 6 S_{10^6} S106​，并指出有效位数
  • 通过本上机题你明白了什么

二、通用程序

def my_fun(n):
    sn, sn1, sn2 = 1/2 * (3/2 - 1/n - 1/(n+1)), 0, 0

    for i in range(2, n+1):
        sn1 += 1/(i**2 - 1)

    for i in range(n, 1, -1):
        sn2 += 1/(i**2 - 1)

    print('The N we use in calculate Sn:', n, '\nThe accurate sum: Sn =', sn,
          '\nFrom large to small sum: Sn1 =', sn1, '\nFrom small to large sum: Sn2 =', sn2)
    print('abs(Sn-Sn1) =', abs(sn-sn1), '\nabs(Sn-Sn2) =', abs(sn-sn2))
    print('-' * 55)


if __name__ == '__main__':
    N = [100, 10000, 1000000]

    for num in N:
        my_fun(num)

三、程序结果

The N we use in calculate Sn: 100 
The accurate sum: Sn = 0.740049504950495 
From large to small sum: Sn1 = 0.7400495049504949 
From small to large sum: Sn2 = 0.7400495049504949
abs(Sn-Sn1) = 1.1102230246251565e-16 
abs(Sn-Sn2) = 1.1102230246251565e-16
-------------------------------------------------------
The N we use in calculate Sn: 10000 
The accurate sum: Sn = 0.7499000049995 
From large to small sum: Sn1 = 0.7499000049995057 
From small to large sum: Sn2 = 0.7499000049995
abs(Sn-Sn1) = 5.662137425588298e-15 
abs(Sn-Sn2) = 0.0
-------------------------------------------------------
The N we use in calculate Sn: 1000000 
The accurate sum: Sn = 0.7499990000005 
From large to small sum: Sn1 = 0.7499990000005217 
From small to large sum: Sn2 = 0.7499990000004999
abs(Sn-Sn1) = 2.1649348980190553e-14 
abs(Sn-Sn2) = 1.1102230246251565e-16
-------------------------------------------------------

四、有效位分析

n==10^2 n==10^4 n==10^6
从小到大 15 ? 15
从大到小 15 13 13

五、结果分析

  • 从有效位可以看出，程序算法对实验的误差具有一定的影响，一个好的算法可以使结果更加精确
  • 以该程序为例，从大到小会出现“大数吃小数”的现象，导致其精度随轮次的增加而相应减少
  • 从小到大顺序可能因为机器问题，导致其误差为 0

第二章

一、问题

  • 给定初值 x 0 x_0 x0​ 及容许误差 ϵ \epsilon ϵ，编制 Newton 法解方程 f ( x ) = 0 f(x) = 0 f(x)=0 根的通用程序
  • 给定方程 f ( x ) = x 3 3 − x = 0 f(x) = \frac{x^3}{3} - x = 0 f(x)=3x3​−x=0，易知其有三个根 x 1 ∗ = − 3 x_1^* = -\sqrt{3} x1∗​=−3 ​， x 2 ∗ = 0 x_2^* = 0 x2∗​=0， x 3 ∗ = 3 x_3^* = \sqrt{3} x3∗​=3 ​
    • 由 Newton 方法的局部收敛性可知存在 δ > 0 \delta \gt 0 δ>0，当 x 0 ∈ ( − δ , δ ) x_0 \in (-\delta, \delta) x0​∈(−δ,δ) 时 Newton 迭代序列收敛于根 x 2 ∗ x_2^* x2∗​，试确定尽可能大的 δ \delta δ
    • 试取若干初始值，观察当 x 0 ∈ ( − ∞ , − 1 ) x_0 \in (-\infty, -1) x0​∈(−∞,−1)， ( − 1 , − δ ) (-1, -\delta) (−1,−δ)， ( − δ , δ ) (-\delta, \delta) (−δ,δ)， ( δ , 1 ) (\delta, 1) (δ,1)， ( 1 , + ∞ ) (1, +\infty) (1,+∞) 时 Newton 序列是否收敛以及收敛于哪一根
  • 通过本上机题，你明白了什么

二、通用程序

from sympy import *


# 牛顿法通用程序
def newton_fun():
    x = Symbol('x')
    epson = eval(input('Enter the absolute error: '))  # 误差
    x0 = eval(input('Enter the initial value: '))      # 初值
    fun = eval(input('Enter the function: '))          # 函数

    derivative = diff(fun, x, 1)                       # 一阶导
    x1 = x0 - fun.subs('x', x0) / derivative.subs('x', x0)

    while abs(x1 - x0) > epson:
        # Newton 迭代格式
        x0, x1 = x1, x1 - fun.subs('x', x1) / derivative.subs('x', x1)
        print(x0, x1)


# 返回收敛值
def discriminant(left, right, x0=0):
    x = Symbol('x')
    epson, fun = 1e-13, x ** 3 / 3 - x
    x0 = (left + right) / 2 if x0 == 0 else x0

    derivative = diff(fun, x, 1)
    x1 = x0 - fun.subs('x', x0) / derivative.subs('x', x0)

    while abs(x1 - x0) > epson:
        x0, x1 = x1, x1 - fun.subs('x', x1) / derivative.subs('x', x1)

    return x1


# 二分法搜索
def search():
    l, r, delta, e = 0.0, 1.0, 0.0, 1e-13
    while abs(delta - (l + r) / 2) > e:
        delta = (l + r) / 2
        res = discriminant(l, r)
        if res == 0:
            l = (l + r) / 2  # 收敛于 0 则将 l 设为中点
        else:
            r = (l + r) / 2  # 不收敛于 0（收敛于根号3） 则将 r 设为中点

    return l  # delta 值可能使函数不收敛于0，故返回 l


if __name__ == '__main__':
    newton_fun()
    print('-' * 50)
    delta = search()
    print('delta =', delta)
    print('-' * 50)

    # 题目所给区间随机选取数值
    arr = [[-100.0, -50.0, -20.0, -10.0, -5.0, -3.0, -1.5], [-0.95, -0.93, -0.89, -0.85, -0.78],
           [-0.77, -0.53, -0.31, 0.27, 0.61], [0.78, 0.81, 0.89, 0.93, 0.97],
           [1.3, 5.1, 11.2, 20.3, 50.2, 100.0]]

    for a in arr:
        for i in a:
            res = discriminant(0, 0, i)
            print('num ', i, ' ==>', res)
        print('-' * 50)

注：sympy 需自行安装文章来源地址https://www.toymoban.com/news/detail-734342.html

三、程序结果

Enter the absolute error: 1e-17
Enter the initial value: 0.5
Enter the function: x**2 - 2*x
-0.250000000000000 -0.0250000000000000
-0.0250000000000000 -0.000304878048780485
-0.000304878048780485 -4.64611473301397e-8
-4.64611473301397e-8 -1.07931905401436e-15
-1.07931905401436e-15 -5.91645678915759e-31
-5.91645678915759e-31 0
--------------------------------------------------
delta = 0.7745966692414186 
 --------------------------------------------------
num  -100.0  ==> -1.73205080756888
num  -50.0  ==> -1.73205080756888
num  -20.0  ==> -1.73205080756888
num  -10.0  ==> -1.73205080756888
num  -5.0  ==> -1.73205080756888
num  -3.0  ==> -1.73205080756888
num  -1.5  ==> -1.73205080756888
--------------------------------------------------
num  -0.95  ==> 1.73205080756888
num  -0.93  ==> 1.73205080756888
num  -0.89  ==> 1.73205080756888
num  -0.85  ==> 1.73205080756888
num  -0.78  ==> -1.73205080756888
--------------------------------------------------
num  -0.77  ==> 0
num  -0.53  ==> 0
num  -0.31  ==> 0
num  0.27  ==> 0
num  0.61  ==> 0
--------------------------------------------------
num  0.78  ==> 1.73205080756888
num  0.81  ==> -1.73205080756888
num  0.89  ==> -1.73205080756888
num  0.93  ==> -1.73205080756888
num  0.97  ==> -1.73205080756888
--------------------------------------------------
num  1.3  ==> 1.73205080756888
num  5.1  ==> 1.73205080756888
num  11.2  ==> 1.73205080756888
num  20.3  ==> 1.73205080756888
num  50.2  ==> 1.73205080756888
num  100.0  ==> 1.73205080756888
--------------------------------------------------

四、结果分析

  • 尽可能大的 δ \delta δ 为 0.7745966692414186
  • ( − ∞ , − 1 ) (-\infty, -1) (−∞,−1) 区间收敛于 x 1 ∗ x_1^* x1∗​， ( 1 , + ∞ ) (1, +\infty) (1,+∞) 区间收敛于 x 3 ∗ x_3^* x3∗​
  • ( − 1 , − δ ) (-1, -\delta) (−1,−δ) 和 ( δ , 1 ) (\delta, 1) (δ,1) 局部收敛于 x 1 ∗ x_1^* x1∗​ 和 x 3 ∗ x_3^* x3∗​
  • ( − δ , δ ) (-\delta, \delta) (−δ,δ) 区间收敛于 x 2 ∗ x_2^* x2∗​
  • 通过本次上机可知，对于多根方程，同一区间上可能局部收敛于不同的根，且迭代序列收敛于某一根有一定区间限制

第三章

一、问题

  • 编制解 n n n 阶线性方程组 A x = b Ax = b Ax=b 的 S O R SOR SOR 方法的通用程序（要求 ∥ x ( k ) − x ( k − 1 ) ∥ ∞ ≤ ϵ \parallel x^{(k)} - x^{(k-1)} \parallel_{\infty} \le \epsilon ∥x(k)−x(k−1)∥∞​≤ϵ）
  • 对于第 39 题中所给线性方程组，取松弛因子 ω i = i 50 ( i = 1 , 2 , ⋯   , 99 ) \omega_i = \frac{i}{50}(i=1, 2,\cdots ,99) ωi​=50i​(i=1,2,⋯,99)，容许误差 ϵ = 1 2 × 1 0 − 5 \epsilon = \frac{1}{2} \times 10^{-5} ϵ=21​×10−5，打印松弛因子，迭代次数，最佳松弛因子及解向量

二、通用程序

import numpy as np


def sor(arr, b, w=0.5, epson=1e-13):
    x = b
    arr_dia = arr.diagonal().reshape((-1, 1))  # 计算矩阵的相应对角矩阵
    # SOR 迭代格式
    x1 = (1 - w) * x + w * (b - arr.dot(x) + arr_dia * x) / arr_dia
    count = 0  # 迭代次数计数器

    e = max(abs(x1 - x))[0]  # ||x1 - x||∞
    while e > epson:
        x, x1 = x1, (1 - w) * x1 + w * (b - arr.dot(x1) + arr_dia * x1) / arr_dia
        e = max(abs(x1 - x))[0]
        count += 1

    return x, count


def sor_test():
    arr = np.array(eval(input('Enter the matrix: ')))
    b = np.array(eval(input('Enter the solution: ')))
    x, count = sor(arr, b)
    print('The solution: ', x)
    print('Iteration counter: ', count)


def sor_test2():
    # 39 题方程组相关矩阵
    arr = [[31, -13, 0, 0, 0, -10, 0, 0, 0],
           [-13, 35, -9, 0, -11, 0, 0, 0, 0],
           [0, -9, 31, -10, 0, 0, 0, 0, 0],
           [0, 0, -10, 79, -30, 0, 0, 0, -9],
           [0, 0, 0, -30, 57, -7, 0, -5, 0],
           [0, 0, 0, 0, -7, 47, -30, 0, 0],
           [0, 0, 0, 0, 0, -30, 41, 0, 0],
           [0, 0, 0, 0, -5, 0, 0, 27, -2],
           [0, 0, 0, -9, 0, 0, 0, -2, 29]]
    arr = np.array(arr)
    b = np.array([[-15], [27], [-23], [0], [-20], [12], [-7], [7], [10]])

    optimal = [0, 1e9, []]  # 最佳松弛因子，迭代次数，解向量
    for i in range(1, 100):
        wi = i / 50
        res, count = sor(arr, b, w=wi, epson=5e-6)

        # 根据迭代次数选取最佳松弛因子及相关解向量
        if optimal[1] > count:
            optimal[0], optimal[1], optimal[2] = wi, count, res

        print('res =', res)
        print('count =', count, '\twi =', wi)
        print('-' * 35)
    print('The optimal wi: ', optimal[0], '\tThe optimal count: ', optimal[1])
    print('The solution: ', optimal[2])


if __name__ == '__main__':
    sor_test2()
    sor_test()

三、程序部分结果

-----------------------------------
res = [[-7.41956679e+306]
 [ 6.76846047e+306]
 [-3.96298666e+306]
 [ 2.72277429e+306]
 [-3.44056652e+306]
 [             inf]
 [            -inf]
 [ 1.01974365e+306]
 [-1.27575693e+306]]
count = 804 	wi = 1.98
-----------------------------------
The optimal wi:  0.94 	The optimal count:  37
The solution:  [[-0.28924816]
 [ 0.34542227]
 [-0.71281942]
 [-0.22061393]
 [-0.43040709]
 [ 0.1542935 ]
 [-0.05783802]
 [ 0.20105187]
 [ 0.29022619]]
-----------------------------------
Enter the matrix: [[1, 2, 1, -2], [2, 5, 3, -2], [-2, -2, 3, 5], [1, 3, 2, 3]]
Enter the solution: [[4], [7], [-1], [0]]
The solution:  [[ 2.]
 [-1.]
 [ 2.]
 [-1.]]
Iteration counter:  1301

四、结果分析

  • 函数 S O R SOR SOR 使用 S O R SOR SOR 迭代格式求解，默认松弛因子为 0.5，默认误差为 1e-13
  • 函数 sor_test 对迭代格式的通用性进行验证，其结果与真实一致
  • 函数 sor_test2 通过多轮迭代求解出最佳松弛因子及解向量
  • 由于得出的解在误差范围内，所以最佳松弛因子由迭代次数最少来确定
  • 当松弛因子大于等于 1.18 时，会出现迭代次数过多以及 inf、nan 等现象？因此，选取松弛因子时应该考虑其大小范围，避免出现错误

到了这里，关于数值分析上机题（上）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处： 如若内容造成侵权/违法违规/事实不符，请点击违法举报进行投诉反馈，一经查实，立即删除！

分享到：
领支付宝红包 赞助服务器费用

相关文章

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包，优惠每天领

二维码1

领取红包

二维码2

领红包