exp1_1.asv

使用matlab软件编制的计算程序

% exp1_1.m --- 算法的数值稳定性实验
% 见 P1 例1 求积分 
%                       I(n) = exp(-1)*int( x^n*exp(x),0,1 )

function try_stable
global n          % 定义全局变量 n ,见后自定义函数 f 中的参量            
N = 20;           % 计算 N 个值

%--------------------------------------------
% [算法1] 用递推公式 
%                        I(k) = 1 - k*I(k-1)
% 取初值                 I0=1-exp(-1)  (约有15位有效数字)     
% [注] 如果取 I0 = 0.6321,结果同书 P2 表1-1,这里初值更精确

I0 = 1-exp(-1);          % 初值
I = zeros(N,1);          % 创建 N x 1 矩阵(即列向量),元素全为零
I(1) = 1-I0;
for k = 2:N
    I(k) = 1 - k*I(k-1);
end
%--------------------------------------------
% [算法2] 用递推公式
%                        I(k-1) = ( 1 - I(k) ) / k
% 取初值                 II(N) = 0 (参见 P5 例3)

II = zeros(N,1);
II(N) = 0;
for k = N:-1:2
    II(k-1) = ( 1 - II(k) ) / k;
end
%--------------------------------------------
%  调用 matlab 高精度数值积分命令 quadl 计算以便比较
%  [注] 该命令以后学习,你现在可模仿使用

III = zeros(N,1);
e_1 = exp(-1);
for k = 1:N
    n = k;                            % 给函数 f 中的参量 n 赋值 
    III(k) = e_1 * quadl(@f,0,1);     % 求函数 f 在[0,1]上的定积分
end
%--------------------------------------------
% 显示计算结果
clc   % 清命令窗
fprintf('\n               算法1结果           算法2结果            精确值')
for k = 1:N,
    fprintf('\nI(%2.0f)  %17.7f   %17.7f   %17.7f',k,I(k),II(k),III(k))
end
% [注]  这里所谓的精确值是指计算显示的数字全是有效数字
%--------------------------------------------
function y = f(x)         % 定义函数
global n                  % 参量 n 为全局变量
y = x.^n.*exp(x);         % ★注意:这里一定要 '点' 运算
return
%--------------------------------------------


% ********   思考题   ********
% 通过上述结果,说明了什么?
% 是不是公式是对的,程序是对的,计算的结果一定是可靠的?
% 造成这种现象的原因是什么? 如何从理论上加以分析?

% ********  你的实验  ********

% 【实验一】  二次方程求根( 见P8 例6 )
% 编写二次方程求根的数值稳定的算法

% ★【实验二】  P11 实验课题(一)
% [说明] 你可模仿该程序来做实验,你只要把下面程序复制为新的 M-文件,
%        去掉开头的 % (参见菜单Text中的 Ctrl_T 与 Ctrl_R命令)
%        再把 '?' 部分改写正确就是一个完整的程序


% myexp1_1.m --- 算法的数值稳定性实验
% 见 P11 实验课题(一) 
%
 function try_stable
 global n a
 N = 20;                  % 计算 N 个值
 a = 0.05;
% %--------------------------------------------
% % [方案I] 用递推公式 
% %         I(k) =  - a*I(k-1) + 1/k
% 
 I0 = log((a+1)/a);                  % 初值
 I = zeros(N,1);          % 创建 N x 1 矩阵(即列向量),元素全为零
 I(1) =-a*I(0)+1;
 for k = 2:N
     I(k) =(-a)*I(k-1)+1/k;
 end
% %--------------------------------------------
% % [方案II] 用递推公式
% %          I(k-1) = ( - I(k) + 1/k ) / a
% 
 II = zeros(N,1);
 if a >= N/(N+1)
     II(N) =(2*a+1)/(2*a(a+1)(N+1));
 else
     II(N) =1/2*(1/((a+1)*(N+1))+1/N);
 end
 for k = N:-1:2
     II(k-1) = 1/a*(-II(k)+1/k);
 end
% %--------------------------------------------
% %  调用 matlab 高精度数值积分命令 quadl 计算以便比较
% 
 III = zeros(N,1);
 for k = 1:N
     n = k;
     III(k) = quadl(@f,0,1);
 end
% %--------------------------------------------
% % 显示计算结果
% 
 clc 
 fprintf('\n               方案I结果           方案II结果            精确值')
 for k = 1:N,
     fprintf('\nI(%2.0f)  %17.7f   %17.7f   %17.7f',k,I(k),II(k),III(k))
 end
% %--------------------------------------------
 function y = f(x)         % 定义函数
 global n a                % 参量 n 为全局变量
 y =x.^n./(a+x);                    % ★注意:这里一定要 '点' 运算
 return
% %--------------------------------------------