第一章
一、问题
设 S N = ∑ j = 2 N 1 j 2 − 1 S_N = \sum_{j=2}^N \frac{1}{j^2 - 1} SN=∑j=2Nj2−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)文章来源:https://www.toymoban.com/news/detail-734342.html
- 编制按从大到小的顺序 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模板网!