实验一算法的数值稳定性实验

% 见 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

% %--------------------------------------------