计算物理学（数值分析）上机实验答案1、误差分析

实验一、误差分析

​ 误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时， 由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法的好坏会影响到数值结果的精度。

一、实验目的

1、通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；

2、通过上机计算，了解误差、绝对误差、误差界、相对误差界的有关概念；

3、通过上机计算，了解舍入误差所引起的数值不稳定性。

二、算法实例

例 1.1

​ 用差商$f^{\prime }(a)\approx \frac{f(a+h)-f(a)}{h}$求 f (x) = ln x 在 x = 3 处导数的近似值。取h= 0.1，h = 0.0001，h =0.000 000 000 000 001 和h=0.000 000 000 000 000 1 分别用 MATLAB 软件计算，取十五位数字计算。

a=3;h=0.1;y=log(a+h)-log(a);yx=y/h
%运行后得yx=0.32789822822991
%将此程序中h改为 0.000 1，运行后得yx=0.33332777790385
%后者比前者好。
%再取h=0.000 000 000 000 001，
%运行后得yx=0.44408920985006
%不如前者好。
%取h=0.000 000 000 000 000 1
%运行后得yx = 0
%算出的结果反而毫无价

例 1.2

​ 分别求方程组 AX = b 在下列情况时的解，其中$A=\left(\begin{array}{cc}1& 1\\ 1& 1.01\end{array}\right)$.
$$\text{ (1) }b=\left(\begin{array}{l}2\\ 2\end{array}\right);\text{ (2) }b=\left(\begin{array}{c}2\\ 2.01\end{array}\right)$$

解： （1） 首先将方程组 AX = b 化为同解方程 X = A⁻¹b ，然后在 MATLAB 工作窗口输入程序

b=[2,2]';A=[1,1;1,1.01]; X=A\b

运行后输出：当$b=\left(\begin{array}{c}2\\ 2.01\end{array}\right)$时， AX = b 的解为$X=\left(\begin{array}{l}2\\ 0\end{array}\right)$ ；

（2）同理可得，当$b=\left(\begin{array}{c}2\\ 2.01\end{array}\right)$时，AX = b 的解为$X=\left(\begin{array}{l}1\\ 1\end{array}\right)$.

例 1.3

​ 计算e的近似值。

解：泰勒级数
$$e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots +\frac{x^n}{n!}+\cdots (-\infty <x<\infty )$$
取 x = 1，得（1.1）式：

$$e=1+1+\frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots +\frac{1}{n!}+\cdots \text{\hspace{1em}(1.1)}$$
这是一个无限过程，计算机无法求到精确值。只能在（1.1）取有限项时计算，再估计误差。如果取有限项

$$s_n(1)=1+1+\frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots +\frac{1}{n!}$$
作为e 的值必然会有误差，根据泰勒余项定理可知其截断误差为
$$e-s_n(1)=\frac{e^{\theta }}{(n+1)!}(0<\theta <1)$$
如果取（1.1）的前九项，输入程序

n=8;s=1;S=1;
for k=1:n;
	s=s*k;
	S=S+1/s;
end
S,R=3/(s*(n+1))
% 运行后结果:
% S = 2.7183	R = 8.2672e-06

因为截断误差为
$$e-s_g(1)=\frac{e^{\theta }}{(8+1)!}<\frac{3}{9!}\approx 8.267196\times 10^{-6}(0<\theta <1)$$
所以 e 的近似值$e^{\approx }{s}_8(1)=1+1+\frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\frac{1}{5!}+\frac{1}{6!}+\frac{1}{7!}+\frac{1}{8!}\approx 2.71828$

例 1.4

​ 取2.718 28作为e的四舍五入近似值时，求其绝对误差和相对误差。

juewu=exp(1)-2.71828	%juewu = 1.8285e-06

例 1.5

​ 计算${\int }_0^{\frac{\pi }{2}}\frac{\sin x}{x}dx$的近似值，并确定其绝对误差和相对误差。

解 因为被积函数$\frac{\sin x}{x}$的原函数不是初等函数，故用泰勒级数求之。

$$\frac{\sin x}{x}=1-\frac{x^2}{3!}+\frac{x^4}{5!}-\frac{x^6}{7!}+\frac{x^8}{9!}+\cdots (-\infty <x<\infty )\text{\hspace{1em}(1.5)}$$
这是一个无限过程， 计算机无法求到精确值。 可用（ 1.5 ）的前四项$1-\frac{x^2}{3!}+\frac{x^4}{5!}-\frac{x^6}{7!}$代替被积函数$\frac{\sin x}{x}$，得

$$y={\int }_0^{\frac{\pi }{2}}\frac{\sin x}{x}dx\approx {\int }_0^{\frac{\pi }{2}}(1-\frac{x^2}{3!}+\frac{x^4}{5!}-\frac{x^6}{7!})dx=\frac{\pi }{2}-\frac{(\frac{\pi }{2}{)}^3}{3\cdot 3!}+\frac{(\frac{\pi }{2}{)}^5}{5\cdot 5!}-\frac{(\frac{\pi }{2}{)}^7}{7\cdot 7!}=\hat{y}$$
根据泰勒余项定理和交错级数收敛性的判别定理，得到绝对误差$R=|y-\hat{y}|<\frac{(\frac{\pi }{2}{)}^9}{9\cdot 9!}=WU$

在 MATLAB 命令窗口输入计算程序如下：

syms x
f=1-x^2/(1*2*3)+x^4/(1*2*3*4*5)-x^6/(1*2*3*4*5*6*7)
y=int(f,x,0,pi/2),y1=double(y)
y11=pi/2-(pi/2)^3/(3*3*2)+(pi/2)^5/(5*5*4*3*2)-(pi/2)^7/(7*7*6*5*4*3*2)
inf=int(sin(x)/x,x,0,pi/2) ,infd=double(inf)
WU =(pi/2)^9/(9*9*8*7*6*5*4*3*2), R =infd-y11

因为运行后输出结果为：y = 1.370 762 168 154 49，yˆ =1.370 744 664 189 38，R = 1.750 396 510 491 47e-005，WU= 1.782 679 830 970 664e-005 < 10−4 . 所以，yˆ 的绝对误差为ε = 10^-4，故$y={\int }_0^{\frac{\pi }{2}}\frac{\sin x}{x}dx\approx 1.3707$ 。 yˆ 的相对误差为${\epsilon }_r=\frac{\epsilon }{|\hat{y}|}=\frac{10^{-4}}{1.3707}<0.0073\%$

例 1.6

​ 设计三种算法求方程2x² +x−15=0在(2, 3) 的一个正根x* 的近似值，并研究每种算法的误差传播情况.

解：设计如下三种算法，然后将计算结果与此方程的精确解 x* = 2.5 比较，观察误差的传播.

算法 1 将已知方程化为同解方程$x=15-2x^2$ 取初值 x₀ = 2 ，按迭代公式$x_{k+1}=15-2x_k^2$依次计算$x_1,x_2,\cdots ,x_n,\cdots$

算法 2 将已知方程化为同解方程$x=\frac{15}{2x+1}$取初值 x₀= 2 ，按迭代公式$x_{k+1}=\frac{15}{2x_k+1}$依次计算$x_1,x_2,\cdots ,x_n,\cdots$

算法 3 将已知方程化为同解方程$x=x-\frac{2x^2+x-15}{4x+1}$取初值 x₀ = 2 ，按迭代公式为$x_{k+1}=x_k-\frac{2x_k^2+x_k-15}{4x_k+1}$依次计算$x_1,x_2,\cdots ,x_n,\cdots$

1. MATLAB 主程序

function [k,juecha,xiangcha,xk]= liti1_6(fun,x0,x1,limax)
% 输入的量--fun是迭代公式，x0是初值, limax是迭代次数和精确值x;
% 输出的量--每次迭代次数k、迭代值xk、绝对误差juecha和相对误差xiangcha，
x(1)=x0;guocheng=[];
for i=1:limax
    x(i+1)=feval(fun,x(i));
    juecha = abs(x(i)-x1);
    xiangcha = juecha /( abs(x(i))+eps);
    xk=x(i);k=i-1;
    guocheng=[guocheng;(i-1),juecha,xiangcha,xk];
end
guocheng

2. MATLAB 调用函数程序及其计算结果

算法 1 fun=@(x)(15-x^2); [k,juecha,xiangcha,xk]= liti1_6(fun,2,2.5,10)

算法 2 fun=@(x)(15/(2*x+1)); [k,juecha,xiangcha,xk]= liti1_6(fun,2,2.5,10)

算法 3 fun=@(x)(x-(2*x^2+x-15)/(4*x+1)); [k,juecha,xiangcha,xk]= liti1_6(fun,2,2.5,10)

​ 运行后将三种算法的前 8 个迭代值$x_1,x_2,\cdots ,x_8$列在一起（见表 1-1），进行比较.将三种算法的$x_1,x_2,\cdots ,x_8$对应的绝对误差和相对误差的值列在一起（见表1-2），进行比较。

表 1-1 例 1.6 中三种算法的计算结果

算 法/迭代次数	算法 1 的迭代结果	算法 2 的迭代结果	算法 3 的迭代结果
> 0	> 2	> 2.000 000 00	> 2.000 000 00
> 1	> 7	> 3.000 000 00	> 2.555 555 56
> 2	> -83	> 2.142 857 14	> 2.500 550 06
> 3	> -13 763	> 2.837 837 84	> 2.500 000 06
> 4	> -378 840 323	> 2.246 963 56	> 2.500 000 00
> 5	> -2.870 4 ×1017	> 2.246 963 56	> 2.500 000 00
> 6	> -1.647 8 ×1035	> 2.321 774 84	> 2.500 000 00
> 7	> -5.430 7 ×1070	> 2.657 901 65	> 2.500 000 00
…	…	…	…
> 99	> -Inf	> 2.500 000 01	> 2.500 000 00

表 1-2 例 1.6 中三种算法计算结果的误差

> 算法	> 算法 1 的误差		> 算法 2 的误差		> 算法 3 的误差
> 迭代次数	> 绝对误差	> 相对误差	> 绝对误差	> 相对误差	> 绝对误差	> 相对误差
> 0	> 0.500 000 00	0.250 000 00	0.500 000 00	> 0.250 000 00	> 0.500 000 00	> 0.250 000 00
> 1	> 4.500 000 00	> 0.642 857 14	0.500 000 00	> 0.166 666 67	> 0.055 555 60	> 0.021 739 13
> 2	> 85.500 000 00	1.030 120 48	0.357 142 86	> 0.1666 666 70	> 0.000 550 10	> 0.000 219 97
> 3	> 13 765.500 00	0.000 100 02	0.337 837 84	> 0.119 047 62	> 0.000 000 06	> 0.000 000 02
> 4	> 378 840 326	1.000 000 01	0.253 036 44	> 0.112 612 62	> 0.000 000 00	> 0.000 000 00
> 5	> 2.870 399 81×1017	> 1	0.230 287 04	> 0.084 345 48	> 0	> 0
> 6	> 1.647 839 01×1035	> 1	0.178 225 16	> 0.076 762 47	> 0	> 0
> 7	> 5.430 746 80×1070	> 1	0.157 901 65	> 0.059 408 39	> 0	> 0
…	…	…	…	…	…	…
> 99	> Inf	> NaN	> 0.000 000 01	> 0.000 000 00	> 0	> 0

例 1.7

​ 求数$x=7^{15}\times (\sqrt{1+8^{-19}}-1)$的近似值。

解 （1）直接用 MATLAB 命令

x=(7^15)*(sqrt(1+8^(-19))-1)	% 运行后输出结果x = 0

问题出现在两个相近的数$\sqrt{1+8^{-19}}$与1相减时，计算机运行程序

sqrt(1+8^(-19))-1	% 运行后输出结果ans = 0

由于计算机硬件只支持有限位机器数的运算，因此在计算中可能引入和传播舍入误差.因为有效数字的严重损失，导致输出$\sqrt{1+8^{-19}}-1$的结果为 0，计算机不能再与数7¹⁵ 继续进行真实的计算，所以，最后输出的结果与 x 的精确值不符。

（2）如果化为$x=7^{15}\times (\sqrt{1+8^{-19}}-1)=\frac{7^{15}\times 8^{-19}}{\sqrt{1+8^{-19}}+1}$再用 MATLAB 命令

x=(7^15)*( (8^(-19))/(sqrt(1+8^(-19))+1))	% 运行后输出结果x = 1.6471e-005

这是因为$\sqrt{1+8^{-19}}-1$化为后$\frac{8^{-19}}{\sqrt{1+8^{-19}}+1}$，计算机运行程序

x= (8^(-19))/(sqrt(1+8^(-19))+1)	% 运行后的结果为x =3.4694e-018

由于有效数字的损失甚少，所以运算的结果3.4694×10^-18 再与7¹⁵ 继续计算，最后输出的结果与 x 的精确值相差无几。

例 1.8

求数$y=\ln (30-\sqrt{30^2-1})$的近似值。

解 （1）直接用 MATLAB 程序

x=30;x1= sqrt(x^2-1)			%运行后输出结果x1 = 29.9833
x=30; x1=29.9833;y=log(x-x1)	%运行后输出结果y = -4.0923

（2）因为$\ln (30-\sqrt{30^2-1})$中的 x=30很大，如果采用倒数变换法
$$z_1=x-\sqrt{x^2-1}=\frac{1}{x+\sqrt{x^2-1}}$$

即
$$\ln (30-\sqrt{30^2-1})=\ln \frac{1}{30+\sqrt{30^2-1}}=-\ln (30+\sqrt{900-1})$$

x=30;y=-log(x+sqrt(x^2-1))	%运行后输出结果y = -4.0941
x=30;y= log(x-sqrt(x^2-1))	%运行后输出结果y = -4.0941

可见，（2）计算的近似值比（1）的误差小。

​ 参加计算的数，有时数量级相差很大.如果不注意采取相应的措施，在它们的加减法运算中，绝对值很小的那个数经常会被绝对值较大的那个数"吃掉"，不能发挥其作用，造成计算结果失真。

例 1.9

请在 16 位十进制数值精度计算机上利用软件 MATLAB 计算下面的两个数x* = 111 111 111 111 111+ 0.1 + 0.3 和 y* = 1 111 111 111 111 111+ 0.1 + 0.3将计算结果与准确值比较，解释计算结果。

解 在 MATLAB 工作窗口输入下面程序

>> x=111111111111111+0.1+0.3, y=1111111111111111+0.1+0.3

运行后输出结果x = 1.111111111111114e+014，y =1.111111111111111e+015

从输出的结果可以看出，x = x* ，而 y ≠ y* .为什么 y * 仅仅比 x* 多一位 1，而 y ≠ y*呢？

​ 这是因为**计算机进行运算时，首先要把参加运算的数写成绝对值小于 1 而"阶码"相同的数，这一过程称为数的"对阶"。**在 16 位十进制数值精度计算机上利用软件 MATLAB 计算这两个数，把运算的数 x* 写成浮点规格化形式为x* = 0.111 111 111 111 111 0 ×1015 + 0.000 000 000 000 000 1×1015 + 0.000 000 000 000 000 3×1015，在 16 位十进制数值精度计算机上，三项的数都表示为小数点后面 16 位数字的数与1015 之积，所以，计算机没有对数进行截断，而是按原来的三个数进行计算。因此，计算的结果 x = x* 。而y* = 0.111 111 111 111 111 1×1016 + 0.000 000 000 000 000 01×1016 + 0.000 000 000 000 000 03×1016三项的数都表示写成绝对值小于 1 而"阶码"都为1016 的数以后，第一项的纯小数的小数点后面有 16 位数字.但是，后两项数的纯小数的小数点后面有 17 位数字，超过了 16 位十进制数值精度计算机的存储量，计算机对后两项的数都进行截断最后一位，即后两项的数都是 16 位机上的零，再进行计算，所以计算结果与实际不符。

三、实验任务

对$n=0,1,2,\cdots ,20$，计算定积分

$$y_n={\int }_0^1\frac{x^n}{x+5}dx$$
算法 1：利用递推公式

$$y_n=\frac{1}{n}-5y_{n-1},n=1,2,\cdots ,20$$
取$y_0={\int }_0^1\frac{1}{x+5}dx=\ln 6-\ln 5\approx 0.182322$

算法 2：利用递推公式
$$y_{n-1}=\frac{1}{5n}-\frac{1}{5}{y}_{n}n=20,19,\cdots ,1$$
注意到

$$\frac{1}{126}=\frac{1}{6}{\int }_0^1{x}^{20}dx\le {\int }_0^1\frac{x^{20}}{x+5}dx\le \frac{1}{5}{\int }_0^1{x}^{20}dx=\frac{1}{105}$$
取$y_{20}\approx \frac{1}{20}(\frac{1}{105}+\frac{1}{126})\approx 0.008730$

y0=0.182322;y20=0.008730;n=20;						% 输入y0,y20是初值, n是迭代次数
y1(1)=y0;y2(n+1)=y20;k=0:n;syms x;
yn=arrayfun(@(n)(int(x^n/(x+5),0,1)),k);			% 准确值计算yn
for i=1:n       y1(i+1)=1/i-5*(y1(i));      end		% 算法一y1
for i=n:-1:1	y2(i)=1/5*(1/i-y2(i+1));	end		% 算法二y2
y=[k',double(yn'),y1',y2']