计算机程序设计艺术习题解答(Excercise 1.2.2-27题)- 求对数标准算法的有限精度误差估计

27. [M25] 此题需要结合正文中计算 \log_{10}x 的标准算法。 用 x_{k}^{'} 表示 x_{k} 基于有限精度的近似值: x(1-\sigma )\leq 10^{n}x_{0}^{'}\leq x(1+\epsilon ) ,用 y_{k} 替代公式(18)中的 (x_{k-1}^{'})^{2},即 (x_{k-1}^{'})^{2}(1-\sigma )\leq y_{k}\leq (x_{k-1}^{'})^{2}(1+\epsilon ) 并且 1\leq y_{k}< 100 ,这里的  \sigma 和 \epsilon 是很小的大于零的常数,表示基于有限精度的近似值的误差上界和下界。若用 \log_{10}^{'}x 表示计算结果,试证明 k 步之后有 \log_{10}x+2\log_{10}{(1-\sigma )}-1/2^{k}< \log_{10}^{'}x\leq \log_{10}x+2\log_{10}(1+\epsilon ) 。

证明:

在此回顾一下正文中计算 \log_{10}x 的标准算法。

\log_{10}x 用二进制表示为 \log_{10}x=n.b_{1}b_{2}b_{3}...=n+b_{1}/2+b_{2}/4+b_{3}/8+...           (17)

假设 x\geq 1 ,若 x < 1, 则首先计算 \log_{10}1/x

x 的整数部分 n 满足 1\leq x/10^{n}<10 。为得到 b_{1}, b_{2}, ... 的值,我们设 x_{0}=x/10^{n} ,对所有的 k\geq 1 ,

若 x_{k-1}^{2}<10 ,则 b_{k}=0,\: x_{k}=x_{k-1}^{2} ;

x_{k-1}^{2}\geq 10 ,则 b_{k}=1,\, x_{k}=x_{k-1}^{2}/10                                                                        (18)

这个过程的有效性来自如下事实:

对 k = 0, 1, 2, ...,有 1\leq x_{k}=x^{2^{k}}/10^{2^{k}(n+b_{1}/2+...+b_{k}/2^{k})}<10 ,                                     (19)

这个事实可以用归纳法证明:当 k = 0 时,显然有 1\leq x_{0}=x/10^{n}<10 ,

k\geq 1 时,若 x_{k-1} 满足条件 1\leq x_{k-1}=x^{2^{k-1}}/10^{2^{k-1}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}<10 ,

此时 x_{k-1}^{2}=x^{2^{k}}/10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})} 满足 1\leq x_{k-1}^{2}<100 ,

当 1\leq x_{k-1}^{2}<10 时,令 b_{k}=0 ,则 x_{k}=x^{2^{k}}/10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1}+0/2^{k})}=x_{k-1}^{2} ,满足 1\leq x_{k}< 10 ;

当 10\leq x_{k-1}^{2}<100 时,令 b_{k}=1,则 x_{k}=x^{2^{k}}/10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1}+1/2^{k})}=x^{2^{k}}/10^{2^{k}(n+b_{1}/2+...b_{k-1}/2^{k-1})+1}=x_{k-1}^{2}/10 ,仍然满足 1\leq x_{k}< 10 。归纳证明完毕。

回到这道习题,我们有 \log_{10}x+2\log_{10}{(1-\sigma )}-1/2^{k}< \log_{10}^{'}x\leq \log_{10}x+2\log_{10}(1+\epsilon ) \Leftrightarrow \: \log_{10}{x(1-\sigma )^{2}/10^{1/2^{k}}}<\log_{10}^{'}x\leq \log_{10}{x(1+\epsilon )^{2}} \Leftrightarrow \; x(1-\sigma )^{2}/10^{1/2^{k}}<x^{'}\leq x(1+\epsilon )^{2}  \Leftrightarrow \: \: x^{2^{k}}(1-\sigma )^{2^{k+1}}/10 < (x^{'})^{2^{k}}\leq x^{2^{k}}(1+\epsilon )^{2^{k+1}}  

由于 \sigma , \epsilon 都是大于零的很小的常数,可以在这里假定 0 <1-\sigma <1, \, \, 1<1+\epsilon <2 , 则 x^{2^{k}}(1-\sigma )^{2^{k+1}}/10 < x^{2^{k}}(1-\sigma )^{2^{k+1}}/(1-\sigma ) = x^{2^{k}}(1-\sigma )^{2^{k+1}-1} ,而 x^{2^{k}}(1+\epsilon )^{2^{k+1}}>x^{2^{k}}(1+\epsilon )^{2^{k+1}}/(1+\epsilon )=x^{2^{k}}(1+\epsilon )^{2^{k+1}-1} ,因此要证明以上的不等式成立,只需证明比它们更加严格的不等式 x^{2^{k}}(1-\sigma )^{2^{k+1}-1}\leq (x^{'})^{2^{k}}\leq x^{2^{k}}(1+\epsilon )^{2^{k+1}-1} 即可。

类似于公式(19),我们有 1\leq x_{k}^{'}=(x^{'})^{2^{k}}/10^{2^{k}(n+b_{1}/2+...+b_{k}/2^{k})}<10 ,即 (x^{'})^{2^{k}}=10^{2^{k}(n+b_{1}/2+...+b_{k}/2^{k})}x_{k}^{'} ,代入上面的式子就得到要证明的不等式为 x^{2^{k}}(1-\sigma )^{2^{k+1}-1}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k}/2^{k})}x_{k}^{'}\leq x^{2^k}(1+\epsilon )^{2^{k+1}-1}​​​​​​​

继续用归纳法。当 k=0 时,根据题面显然有 x(1-\sigma )\leq 10^{n}x_{0}^{'}\leq x(1+\epsilon ) ,

k\geq 1 时,若 x_{k-1}^{'} 满足 x^{2^{k-1}}(1-\sigma )^{2^{k}-1}\leq 10^{2^{k-1}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}x_{k-1}^{'}\leq x^{2^{k-1}}(1+\epsilon )^{2^{k}-1} ,对上式各个部分分别取平方,得到 x^{2^{k}}(1-\sigma )^{2^{k+1}-2}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}(x_{k-1}^{'})^{2}\leq x^{2^{k}}(1+\epsilon )^{2^{k+1}-2}

先看左边的不等式 x^{2^{k}}(1-\sigma )^{2^{k+1}-2}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}(x_{k-1}^{'})^{2} ,两边同时乘以 (1-\sigma ) ,且根据题面有 (x_{k-1}^{'})^{2}(1-\sigma )\leq y_{k} ,\Rightarrow x^{2^{k}}(1-\sigma )^{2^{k+1}-1}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}(x_{k-1}^{'})^{2}(1-\sigma )\leq 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}y_{k}

由于 1\leq y_{k}< 100 ,当 y_{k}< 10 时,令 b_{k}=0,\, \, x_{k}^{'}=y_{k} ,可得 x^{2^{k}}(1-\sigma )^{2^{k+1}-1}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1}+0/2^{k})}y_{k}= 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1}+b_{k}/2^{k})}x_{k}^{'} ;当 10\leq y_{k}< 100 时,令 b_{k}=1, \, \, x_{k}^{'}=y_{k}/10 ,即 y_{k}=10x_{k}^{'} ,可得 x^{2^{k}}(1-\sigma )^{2^{k+1}-1}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}y_{k}= 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1}+1/2^{k})}x_{k}^{'}= 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1}+b_{k}/2^{k})}x_{k}^{'}

。综合起来,可得 x^{2^{k}}(1-\sigma )^{2^{k+1}-1}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k}/2^{k})}x_{k}^{'} 。

再看右边的不等式 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}(x_{k-1}^{'})^{2}\leq x^{2^{k}}(1+\epsilon )^{2^{k+1}-2} ,两边同时乘以 (1+\epsilon )

且根据题面有 y_{k}\leq (x_{k-1}^{'})^{2}(1+\epsilon ) ,可以推得10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}y_{k}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k-1}/2^{k-1})}(x_{k-1}^{'})^{2}(1+\epsilon )\leq x^{2^{k}}(1+\epsilon )^{2^{k+1}-1}

和上面的推导类似,可得 10^{2^{k}(n+b_{1}/2+...+b_{k}/2^{k})}x_{k}^{'}\leq x^{2^{k}}(1+\epsilon )^{2^{k+1}-1} 。结合起来,就得到了要证明的不等式 x^{2^{k}}(1-\sigma )^{2^{k+1}-1}\leq 10^{2^{k}(n+b_{1}/2+...+b_{k}/2^{k})}x_{k}^{'}\leq x^{2^k}(1+\epsilon )^{2^{k+1}-1}

 。 反推回去,就证明了 \log_{10}x+2\log_{10}{(1-\sigma )}-1/2^{k}< \log_{10}^{'}x\leq \log_{10}x+2\log_{10}(1+\epsilon ) 。这个式子也说明,随着 k 的增大,\log_{10}^{'}x 越来越接近于 \log_{10}x​​​​​​​ 。

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
To plot C versus S for 0 ≤ S ≤ 200, we can use the following MatLab code: ```matlab E = 100; % Exercise price T = 1; % Expiry date r = 0.05; % Risk-free interest rate sigma = 0.3; % Volatility S = linspace(0, 200, 1000); C = zeros(size(S)); for i = 1:numel(S) C(i) = EuropeanCall(S(i), E, r, sigma, T); end plot(S, C); xlabel('S'); ylabel('C'); title('European Call Option Price'); ``` This code first defines the parameters of the option, and then generates a range of values for S between 0 and 200 using the `linspace` function. For each value of S, the code calculates the option price using the `EuropeanCall` function (which can be obtained from various sources). Finally, the code plots the option price as a function of S using the `plot` function. To plot the value of C(S, t) at different values of t between t = 0 and t = T, we can modify the above code as follows: ```matlab E = 100; % Exercise price T = 1; % Expiry date r = 0.05; % Risk-free interest rate sigma = 0.3; % Volatility S = linspace(0, 200, 1000); t = linspace(0, T, 1000); C = zeros(numel(S), numel(t)); for i = 1:numel(S) for j = 1:numel(t) C(i, j) = EuropeanCall(S(i), E, r, sigma, T - t(j)); end end surf(S, t, C); xlabel('S'); ylabel('t'); zlabel('C'); title('European Call Option Price'); ``` This code generates a grid of values for S and t using `linspace`, and then calculates the option price for each combination of S and t using a nested loop. The option price is stored in a matrix `C`, and is plotted as a surface using the `surf` function. As t approaches T, the option price converges to the intrinsic value of the option, which is max(S - E, 0) for a call option. This is because as the expiry date approaches, the option has less time to move in the money, and its time value decreases.

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值