数值计算方法
数值计算方法知识面涉及微积分，线性代数；运用编程的方法来解决关于数值计算的问题，其中重点讨论如何最小化误差，一些方程的数值解法，以及插值和拟合问题。其中的知识可以作为数据挖掘的基础。
原文链接 https://blog.gongyan.me/2017/04/shu-zhi-ji-suan-fang-fa/ 转载请注明，有问题欢迎联系我 gongyanc@outlook.com
习题0
执行>> 5.1-5-0.1和>> 1.5-1-0.5,给出执行结果，并简要分析一下产生现象的原因
代码：
>>x1=5.1-5-0.1
>>x1=
>>-3.608224830031759e-16
>>x2=1.5-1-0.5
>>x2=
>>0
原因：
浮点数表示时数字时位数有限，以双精度浮点数为例，共64位。
对于5.1，符号位一位，指数位11，整数部分占3位，小数部分0.1不能精确表示只能近似用49位表示。5.1-5-0.1，5被精确表示，0.1则用52位近似表示。故5.1中的小数部分0.1在计算机中比0.1小，最后结果为负。而x2=1.5-1-0.5中，所有的数都被精确表示了，故其结果为0.
习题1
设$I_n = \int_0^1 \frac{x^n}{5+x} dx$
(1) 从$I_0$尽可能精确的近似值出发，利用递推式$I_n = -5I_{n-1} + \frac{1}{n}\ \ \ \ (n = 1, 2, 3, ... ,20)$,计算$I_{20}$的近似值；
(2) 从$I_{20}$较粗糙的估计值出发，利用递推式$I_{n-1} = -\frac{1}{5}I_n + \frac{1}{5n}\ \ \ \ (n = 20, 19, ... , 1)$计算$I_0$的近似值；
(3) 分析结果的可靠性以及出现这种现象的原因.
(1)
(一) 程序设计
1 算法
$$I_n = -5I_{n-1} + \frac{1}{n} \ \ \ \ (n = 1, 2, 3, ... , 20)$$
2 说明事项
- $I_0$必须取得精确值，利用$I_n = \int_0^1 \frac{x^n}{5+x} dx$计算$I_0$时，会得到ln6-ln5，此处两数相近不宜相减，应换成ln(6/5)
- 计算出$I_0$后利用循环即可计算出$I_{20}$
(二) 源程序
MATLAB:
function [ y ] = test1_1( n )
y=log(1.2);
for i=1:1:n
y=-5*y+1/i;
end
Python:
import math
def test1_1(n):
y = math.log(1.2, math.e)
for i in range(1, n+1, 1):
y = -5 * y + 1 / i
print("%e" % y)
if __name__ == "__main__":
test1_1(20)
(三) 实验结果
>> test1_1(20)
>> ans =
4.242637044921560e-03
(2)
(一) 程序设计
1 算法
$$I_{n-1} = -\frac{1}{5}I_n + \frac{1}{5n} \ \ \ \ (n = 20, 19, ... , 1)$$
2 说明事项
- 将上一题中的结果去除尾数变为4.0e-03，得到$I_{20}$粗糙估计值
- 计算出$I_{20}$后利用循环即可计算出$I_0$
(二) 源程序
function [ y ] = test1_2
y=4.0e-03;
for i=20:-1:1;
y=y/(-5)+1/(5*i);
end
(三) 实验结果
>> test1_2
>> ans =
1.823215567939546e-01
与log(1.2)的值相同
(3)
对于第一个公式
$$I_n = -5I{n-1} + \frac{1}{n}\ \ \ \ (n = 1, 2, 3, ... , 20)$$
分析其误差传递
$$\varepsilon(I^_n) = |I_n - I^n| = -5|I - I^_{n-1}| = (-5)^{20} \varepsilon(I^0)$$
20次之后其误差惊人，其结果不可取
对于第二个公式
$$I = -\frac{1}{5}I_n + \frac{1}{5n} \ \ \ \ (n = 20, 19, ... , 1)$$
分析其误差传递
$$\varepsilon(I^{n-1}) = |I - I^_{n-1}| = -\frac{1}{5}|I_n - I^_n| = (-\frac{1}{5})^{20} \varepsilon(I^_N)$$
20次之后误差很小
习题2
用下列方法求方程$e^x + 10x - 2 = 0$的近似根，要求误差不超过$\frac{1}{2} \times 10^-3$,并比较计算量
(1) 在区间[0, 1]上用二分法；
(2) 取初值$x_0 = 0$并用迭代过程$x_{k+1} = \frac{2-e^{x_k}}{10}\ \ \ \ (k = 0, 1, 2, ...)$；
(3) 取初值$x_0 = 0$用牛顿法.
(1)
(一) 程序设计
1 算法
二分法通过缩短短区间，一步步逼近方程的根
2 说明事项
此代码只适用于本程序
(二) 源程序
function [ x ] = test2_1
a=0;
b=1;
i=1;
while 1
x=(a+b)/2;
fa=exp(a)+10*a-2;
fb=exp(b)+10*b-2;
fx=exp(x)+10*x-2;
if (abs(fa-fb)<1e-3);
break;
end
if (fa*fx>0)
a=x;
else
b=x;
end
i++;
end
i
(三) 试验结果
>>test2_1
>>i=1.50000000000000e+001
>>ans=9.05456542968750e-002
(2)
(一) 程序设计
1 算法
迭代公式$x_{k+1} = \frac{2-e^{x_k}}{10}\ \ \ \ (k = 0, 1, 2, ...)$
(二) 源程序
function [ y ] = fun2_1( x )
y=x;
i=1;
while 1
y1=y;
y=0.1*(2-exp(y1));
if (abs(y-y1)<5e-4)
break;
end
i++;
end
i
(三) 试验结果
>>test2_2(0)
>>i=4.00000000000000e+000
>>ans=9.05126166743651e-002
(3)
(一) 程序设计
1 算法
牛顿迭代公式$x_{k+1} = x_k - \frac{x_k}{f'(x_k)}$
2 说明事项
- 对于定义的函数f(x), 为了方便计算，我们用数值微分来近似计算一阶导数的值：
$$f'(x_k) \approx \frac{f(x_k+h) - f(x_k)}{h}$$
其中h是预先取定的小步长； - 因为计算公式中有除法，必须对除数为0的情况作异常处理.
(二) 源程序
function [ y ] = fun2_3( x )
y=x;
i=1;
while 1
y1=y;
y=y1-(exp(y1)+10*y1-2)/(exp(y1)+10);
if (abs(y-y1)<5e-4)
break;
end
i++;
end
i
(三) 试验结果
>>test2_3(0)
>>i=2.00000000000000e+000 %i为迭代次数
>>ans=9.05251085833896e-002 %ans为根
