自适应求积算法 MatLab版

%自适应simpson算法
clc
clear
close all;
syms x
f(x)=100/(x^2)*sin(10/x);%目标函数
a=1;b=3;%初始区间
TOL=1E-4;%误差容限
[s,fcount]=acsimpson1(f,a,b,TOL);

t = 1:0.001:3;
plot(t,f(t))
hold on
disp("自适应求积结果:")
result = vpa(s);
fprintf("result= %f",result);
% 画出f(x)图像
function [s,fcount]=acsimpson1(f,a,b,tol)
%此方法为自适应simpson方法求解定积分
%f为函数，a，b分别为积分区间，tol为精度
%s为积分值，fcount为计算函数值次数

%一下为准备输入数据
f = fcnchk(f,'vectorized');       %定义内联函数
h=b-a;             %步长
qu=[a a+h/2 b];    %三个节点
% qu(qu==0)=2^-1024; %避免计算函数值出现NAN的情况
y=feval(f,qu);     %计算函数节点处的值
s=y*[1 4 1]'*h/6;  %Simpson积分
fcount=0;          %定义初始计数器
tol=tol*10;        %调整精度

[s,fcount]=subsim(f,s,qu,y,fcount,tol); %递归函数
fcount=fcount+3;   %计数器加上准备输入数据的三次计算
e