exp1_1.m

使用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 = ?;
% %--------------------------------------------
% % [方案I] 用递推公式 
% %         I(k) =  - a*I(k-1) + 1/k
% 
% I0 = ?;                  % 初值
% I = zeros(N,1);          % 创建 N x 1 矩阵(即列向量),元素全为零
% I(1) = ?;
% for k = 2:N
%     I(k) = ?;
% end
% %--------------------------------------------
% % [方案II] 用递推公式
% %          I(k-1) = ( - I(k) + 1/k ) / a
% 
% II = zeros(N,1);
% if a >= N/(N+1)
%     II(N) = ?;
% else
%     II(N) = ?;
% end
% for k = N:-1:2
%     II(k-1) =  ?;
% 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 = ?;                    % ★注意:这里一定要 '点' 运算
% return
% %--------------------------------------------