首先自定义了一个算符o_plus：当a+b属于[-1,1]时，o_plus(a,b)=a+b；若a+b大于1时，o_plus(a,b)=1；若a+b小于-1时，o_plus(a,b)=-1。例如y=o_plus(x,0)的图像如图1所示。
二重积分的被积函数intexp(x,y)=phi(y)*p1(x)*p2(y),其中只有第一项phi(y)的表达式中引入了自定义算符o_plus,phi(y)函数图像如图2，曲面intexp(x,y)的图像如图3-1和3-2所示，由于引入了特殊算符，曲面在整个被积区间上并非是光滑的（导数不连续），经过分析，此时切面y=-17.4195将该曲面分成了两部分，而这两部分都是光滑的，这可以从phi(y)的图形(图2)以及图3-2中看出.
然后我使用integral函数用先积x后y的积分顺序、先y后x的积分顺序以及分段积分求和3种方法计算，积分区间如图4阴影部分所示，区间受直线x=30cos(b)控制。在不同的b值下获得了三种方法下的积分结果，如图5.
在我看来，分段积分的结果应该是正确的，而在不分段的情况下两种积分顺序的结果都与分段的结果偏离，所以针对这样的非光滑曲面能否用不分段的方法直接积分出较为精确的结果？因为针对图6这样更复杂的情形分段的方法实在过于复杂。
我也试过quadgk函数(高斯积分)，其结果与integral函数一样，我的问题是：针对非光滑曲面（二元函数）能否用不分段的方法直接得到二重积分较为精确的结果，应该怎样操作？附我的代码。
clc,clear;close;
Hc=1;
Hm=30;
f=100000;
k=80000;
cf=2*pi*f/(k*Hc);
syms x y a b
uy1=acos(y/Hm);
uy2=2*pi-acos(y/Hm);
phi=o_plus(1/cf*int(Hm*cos(a)/Hc-y,a,uy1,uy2),1);%phi为y的函数，被积函数的第一项
subplot(1,2,1)
fplot(phi,[-Hm,Hm],'-','LineWidth',2)%画出phi(y)的图像
xlabel('y')
ylabel('phi(y)')
y_p=double(vpasolve(1/cf*int(Hm*cos(a)/Hc-y,a,uy1,uy2)+1==-1,y));%求出分段的位置
disp(['y_p的位置是:',num2str(y_p)])
intexp=phi*px*py;%被积函数表达式
intexp_mf=matlabFunction(intexp);%将被积函数转换为函数句柄的形式
A_yx=@(b) integral(@(y) integral(@(x) intexp_mf(x,y),Hm*cos(b),-y),-Hm,-Hm*cos(b),'ArrayValued',true);%积分顺序为先y后x，积分区间为b的函数
A_xy=@(b) integral(@(x) integral(@(y) intexp_mf(x,y),-Hm,-x),Hm*cos(b),Hm,'ArrayValued',true);%积分顺序为先x后y，积分区间为b的函数
a1=matlabFunction((1/cf*int(Hm*cos(a)/Hc-y,a,uy1,uy2)+1)*px*py);%分段积分，第一部分的被积函数
a2=matlabFunction(-px*py);%分段积分，第二部分的被积函数
A_fd_1=@(b) integral(@(y) integral(@(x) a1(x,y),Hm*cos(b),-y),-Hm,y_p,'ArrayValued',true);%分段积分，第一部分，积分区间为b的函数
A_fd_2=@(b) integral(@(y) integral(@(x) a2(x,y),Hm*cos(b),-y),y_p,-Hm*cos(b),'ArrayValued',true);%分段积分，第二部分，积分区间为b的函数
bb=1.5*pi:0.05*pi:2*pi;%b的取值，决定积分区域
BB=zeros(length(bb));
for i=1:length(bb)
BB(1,i)=A_yx(bb(i));%先积y后x的积分方法
BB(2,i)=A_xy(bb(i));%先积x后y的积分方法
BB(3,i)=A_fd_1(bb(i))+A_fd_2(bb(i));%分段积分求和的方法
end
subplot(1,2,2)
plot(bb,BB(1,:),'-','LineWidth',2)
hold on
plot(bb,BB(2,:),'-','LineWidth',2)
hold on
plot(bb,BB(3,:),'-','LineWidth',2)
xlabel('b')
ylabel('积分结果')
legend('先积x后y','先积y后x','分段积分')
%%%引用的函数%%%
function output=o_plus(a,b)%引入的特殊算符
output=min(max(a + b, -1), 1);
end
function output=px(~)%x分布的概率密度函数
syms x
sigma_x=20;%x分布的标准差
a=20;%x分布的平均值
output=1/((2*pi)^0.5*sigma_x)*exp(-(x-a)^2/(2*(sigma_x)^2));
end
function output=py(~)%y分布的概率密度函数
syms y
sigma_y=20;%y分布的标准差
a=-20;%y分布的平均值，值为x的相反数
output=1/((2*pi)^0.5*sigma_y)*exp(-(y-a)^2/(2*(sigma_y)^2));
end
