用编程实现计算方法/数值分析的概念知识（一）

计算方法/数值分析实验代码

计算方法是个神奇的课程，别的课程都是用计算机做运算，这个课是在教计算机做运算，不得不佩服前人的智慧，能把求导，积分这些玄之又玄的运算找到循环的解决方法。敲完这些代码感觉收获还是挺大的。

对时间复杂度有了新的认识。以前算开方直接就调用sqrt()函数了，从来没考虑过他的实现方法和函数本身的时间复杂度。
计算机的存在就是为了把人从繁杂的计算里解放出来，这些底层的运算实现带给我的震撼比那些精巧的数据结构还要大。

数值计算的误差

对$n=0,1,...,20$，采用下面两种递推公式计算定积分:
$$ I_n=\int^1_0 \frac{x^n}{x+5}dx $$
并根据计算结果分析这两种算法的稳定性。
算法一：利用递推公式：
$$ I_n=\frac{1}{n}-5I_{n-1}(n=1,2,...,20) $$
取
$$ I_0 = \int_0^1 \frac{1}{x+5}dx = ln6-ln5 \approx 0.182322 $$
算法二：利用递推公式：
$$ I_{n-1}=\frac{1}{5}( \frac{1}{n} -I_n ) (n=20,19,...,1) $$
注意到：
$$ \frac{1}{126}=\frac{1}{6}\int^1_0 x^{20} \leqslant \int_0^1 \frac {x^{20}}{x+5}dx \leqslant\frac{1}{5} \int^1_0 x^{20}dx = \frac{1}{105} $$
取
$$ I_{20} \approx \frac{1}{2}(\frac{1}{105}+\frac{1}{126}) \approx 0.008730 $$

解：
算法一应该是经典的蝴蝶效应，$n=0$时只有$10^{-6}$数量级的截断误差，迭代20次以后，其误差完全失控了。

from math import log
n = 0               # 迭代次数

# 两种赋值方式
# list_ = [log(6)-log(5)]         # 保留小数：最后计算结果非常精确
list_ = [0.182322]         # 不保留小数，最后计算结果误差极大

while(n<=20):
    print(n,"\t",list_[n])
    n+=1
    l_n = (1/n)-5*list_[n-1]
    list_.append(l_n)

# print(list_[20])

而算法二看似误差取的非常的激进，但经过20次迭代衰减后，结果跟真实结果相差无几。

n = 20              
list_=[0.008730]
while n>0:
    t = 20-n        # 迭代次数
    print(n,"\t",list_[t])
    n-=1
    l_n=(1/5)*(1/(n+1)-list_[t])
    list_.append(l_n)
    
    

print(n,"\t",list_[20])

非线性方程求根

用牛顿迭代法求方程$f(x)=x^3-x-1=0$在$x_0=1.5$附近的根，要求精确到$10^-4$，输出每次的迭代结果并统计所用的迭代次数。
解：根据牛顿迭代法的运算步骤，进行分析：
原方程为$x^3-x-1$，则其一阶导数为$3x^2-1$，带入公式即可。
```
from sys import maxsize

x_ = [1.5]      # 设初始值为1.5
n = 0           # 迭代次数
err = maxsize
while err>=0.00005:
 print(n, "\t", x_[n])
 n+=1
 x_n = x_[n-1]-(x_[n-1]**3-x_[n-1]-1)/(3*(x_[n-1]**2)-1)
 x_.append(x_n)
 err = abs(x_[n]-x_[n-1])

print(n, "\t", x_[n])
```
粘一下自己的结果输出：
n $x_n$
0 1.5
1 1.3478260869565217
2 1.325200398950907
3 1.3247181739990537
4 1.3247179572447898
在第4次时收敛。

n	$x_n$
0	1.5
1	1.3478260869565217
2	1.325200398950907
3	1.3247181739990537
4	1.3247179572447898

用埃特金加速法求第一题中的方程的根，并与第一题的结果比较，哪种方法的收敛速度更快。
解：埃特金加速法的求解方法跟牛顿迭代法相似，只是多了一步运算步骤。
对变量符号做如下定义：

x_[n]: 为x_n的根
x2_[n]: 为x_n近似根
x3_[n]: 为x_n校正根

from sys import maxsize


n=0     # 迭代次数
x_=[1.5]     
# 初始值为1.5
x2_=[-1]
x3_=[-1]
err = maxsize # 用于表示误差
print('n', "\t", 'x2_[n]', "\t", 'x3_[n]', "\t", 'x_[n]')

while err>=0.00005:
 print(n, "\t", x2_[n], "\t", x3_[n], "\t", x_[n])
 n+=1
 x2_n = x_[n-1]**3-1
 x3_n = x2_n**3-1
 x_n  = x3_n-((x3_n-x2_n)**2)/(x3_n-2*x2_n+x_[n-1])
 x_.append(x_n)
 x2_.append(x2_n)
 x3_.append(x3_n)
 err = abs(x_[n]-x_[n-1])

print(n, "\t", x2_[n], "\t", x3_[n], "\t", x_[n])