西安电子科技大学数值分析第四章实验题第1题

第一题

复化梯形代码

clc; clear;

% 积分函数

syms x

eps = 1e-20;

f = @(x) sqrt(abs(x)) * log(abs(x) + eps);


% 积分区间和精度

a = 0;

b = 1;

e = 0.5*1e-3;


% 初始化变量

err = 1;

estimates = [];

width = b - a;

points = [a, b];

new_pts = [];

sum_fx = 0;


% 主循环

digits(8)

format long

iter = 0;

est = vpa(width / 2 * (f(a) + f(b)));

estimates = [estimates, est];

fprintf('K=%d, 2^K=%d, TK=%0.8f\n', iter, 2^iter, est)


while(err > e)

% 计算新的积分点和 sum_fx

for i = 1:length(points)-1

new_pt = (points(i) + points(i+1)) / 2;

sum_fx = sum_fx + f(new_pt);

new_pts = [new_pts, new_pt];

end

% 更新已选的积分点

points = sort([points, new_pts]);

% 更新迭代次数

iter = iter + 1;

% 更新积分估计值

est = vpa(0.5 * (est + width * sum_fx));

estimates = [estimates, est];

% 更新误差

err = abs(est - estimates(end-1));

fprintf('K=%d, 2^K=%d, TK=%0.8f, epsilon =%0.8f\n', iter, 2^iter, est, err)

% 更新值

width = width / 2;

sum_fx = 0;

new_pts = [];

end


final_value = est;

fprintf('复化梯形求积公式的结果为：%0.8f\n', final_value)

运行结果

K=0, 2^K=1, TK=0.00000000
K=1, 2^K=2, TK=-0.24506454, epsilon =0.24506454
K=2, 2^K=4, TK=-0.35810406, epsilon =0.11303952
K=3, 2^K=8, TK=-0.40809004, epsilon =0.04998598
K=4, 2^K=16, TK=-0.42947458, epsilon =0.02138455
K=5, 2^K=32, TK=-0.43838949, epsilon =0.00891490
K=6, 2^K=64, TK=-0.44203068, epsilon =0.00364120
K=7, 2^K=128, TK=-0.44349365, epsilon =0.00146297
K=8, 2^K=256, TK=-0.44407364, epsilon =0.00057998
K=9, 2^K=512, TK=-0.44430104, epsilon =0.00022740
复化梯形求积公式的结果为：-0.44430104

复化辛普森代码

clc; clear;


% 积分函数

syms x

eps = 1e-20;

f = @(x) sqrt(abs(x)) * log(abs(x) + eps);


% 积分区间和精度

a = 0;

b = 1;

e = 0.5 * 1e-3;


% 初始化变量

err = 1;

estimates = [];

width = b - a;

points = [a, b];

n = 1; % 初始分段数

Sn = ComplexSimpson(a, b, n);


% 主循环

digits(8)

format long

iter = 0;

estimates = [estimates, Sn];

fprintf('K=%d, 2^K=%d, TK=%0.8f\n', iter, 2^iter, Sn)


while(err > e)

% 计算新的积分点和复化辛普森积分

new_pts = linspace(a, b, 2^n+1);

Sn = ComplexSimpson(a, b, 2^n);

% 更新迭代次数

iter = iter + 1;

% 更新误差

err = abs(Sn - estimates(end));

fprintf('K=%d, 2^K=%d, TK=%0.8f, epsilon=%0.8f\n', iter, 2^iter, Sn, err)

% 更新值

points = new_pts;

n = n * 2;

estimates = [estimates, Sn];

end


final_value = Sn;

fprintf('复化辛普森求积公式的结果为：%0.8f\n', final_value)



function result = fun(x)

eps = 1e-10;

result = sqrt(abs(x)) * log(abs(x) + eps);

end


function Sn = ComplexSimpson(xL, xR, n)

% xL, xR：积分区间左右端点，n：分段数

% Sn ：复化辛普森积分结果

h = (xR - xL) / n; % 求步长

Sn = 0;

x = xL:h:xR;

for k = 2:n+1

Sn = Sn + 4 * fun(x(k) - 0.5 * h) + 2 * fun(x(k));

end

Sn = Sn + fun(xL) - fun(xR);

Sn = vpa(Sn * h / 6, 8);

end

运行结果

K=0, 2^K=1, TK=-0.32675271
K=1, 2^K=2, TK=-0.39578390, epsilon=0.06903119
K=2, 2^K=4, TK=-0.42475203, epsilon=0.02896813
K=3, 2^K=8, TK=-0.44136112, epsilon=0.01660909
K=4, 2^K=16, TK=-0.44437684, epsilon=0.00301572
K=5, 2^K=32, TK=-0.44444442, epsilon=0.00006758
复化辛普森求积公式的结果为：-0.44444442