数值计算方法知识面涉及微积分，线性代数；运用编程的方法来解决关于数值计算的问题，其中重点讨论如何最小化误差，一些方程的数值解法，以及插值和拟合问题。

习题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小，最后结果为负。而x2x_2=1.5-1-0.5中，所有的数都被精确表示了，故其结果为0.

习题1

设In=∫01xn5+xdxI_{n}=\int_{0}^{1} \frac{x^{n}}{5+x} dx
(1) 从I0I_0尽可能精确的近似值出发，利用递推式In=−5In−1+1n(n=1,2,3,...,20)I_n = -5I_{n-1} + \frac{1}{n} \quad (n = 1, 2, 3, ... ,20),计算I20I_{20}的近似值；
(2) 从I20I_{20}较粗糙的估计值出发，利用递推式In−1=−15In+15n(n=20,19,...,1)I_{n-1} = -\frac{1}{5}I_n + \frac{1}{5n} \quad (n = 20, 19, ... ,1)计算I0I_0的近似值；
(3) 分析结果的可靠性以及出现这种现象的原因.

问题(1)

程序设计

1 算法

In=−5In−1+1n(n=1,2,3,...,20)I_n = -5I_{n-1} + \frac{1}{n} \quad (n = 1, 2, 3, ... , 20)

2 说明事项

  • I0I_0必须取得精确值，利用In=∫01xn5+xdxI_n = \int_0^1 \frac{x^n}{5+x} dx计算I0I_0时，会得到ln6−ln5ln6-ln5，此处两数相近不宜相减，应换成ln(6/5)ln(6/5)
  • 计算出I0I_0后利用循环即可计算出I20I_{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 算法

In−1=−15In+15n(n=20,19,...,1)I_{n-1} = -\frac{1}{5}I_n + \frac{1}{5n} \quad (n = 20, 19, ... , 1)

2 说明事项

  • 将上一题中的结果去除尾数变为4.0e-03，得到I20I_{20}粗糙估计值
  • 计算出I20I_{20}后利用循环即可计算出I0I_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)

对于第一个公式

In=−5In−1+1n(n=1,2,3,...,20)I_n = -5I{n-1} + \frac{1}{n}\quad (n = 1, 2, 3, ... , 20)

分析其误差传递

ε(In∗)=∣In−In∗∣=−5∣In−1−In−1∗∣=(−5)20ε(I0∗)\varepsilon(I^*_n) = |I_n - I^*_n| = -5|I_{n-1} - I^*_{n-1}| = (-5)^{20} \varepsilon(I^*_0)

20次之后其误差惊人，其结果不可取.
对于第二个公式

In−1=−15In+15n(n=20,19,...,1)I_{n-1} = -\frac{1}{5}I_n + \frac{1}{5n} \quad (n = 20, 19, ... , 1)

分析其误差传递

ε(In−1∗)=∣In−1−In−1∗∣=−15∣In−In∗∣=(−15)20ε(IN∗)\varepsilon(I^*_{n-1}) = |I_{n-1} - I^*_{n-1}| = -\frac{1}{5}|I_n - I^*_n| = (-\frac{1}{5})^{20} \varepsilon(I^*_N)

20次之后误差很小

习题2

用下列方法求方程ex+10x−2=0e^x + 10x - 2 = 0的近似根，要求误差不超过12×10−3\frac{1}{2} \times 10^{-3},并比较计算量
(1) 在区间[0, 1]上用二分法；
(2) 取初值x0=0x_0 = 0并用迭代过程xk+1=2−exk10(k=0,1,2,...)x_{k+1} = \frac{2-e^{x_k}}{10}\quad (k = 0, 1, 2, ...)；
(3) 取初值x0=0x_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 算法

xk+1=2−exk10(k=0,1,2,...)x_{k+1} = \frac{2-e^{x_k}}{10} \quad (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 算法

xk+1=xk−xkf′(xk)x_{k+1} = x_k - \frac{x_k}{f'(x_k)}

2 说明事项

  • 对于定义的函数f(x)，为了方便计算，我们用数值微分来近似计算一阶导数的值： f′(xk)≈f(xk+h)−f(xk)hf'(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为根