用matlab求定积分的上限b的值定积分l=∫(x^2+y^2+z^2)^(1/2)*xd(t)上限为b下限为0;其中z=-R2*sin(t);y=R2*cos(t);x=R2*cos(t)/a+(a^2+1)^(1/2)*cos(t)*sin(t)/(a*(R1^2-R2^2*cos(t)^2);已知l=10；R2=5;R1=10;a=1；积分下限为0；求上

Question

用matlab 求定积分的上限b的值
定积分l=∫(x^2+y^2+z^2)^(1/2)*xd(t)上限为b下限为0;其中z=-R2*sin(t);y=R2*cos(t);x=R2*cos(t)/a+(a^2+1)^(1/2)*cos(t)*sin(t)/(a*(R1^2-R2^2*cos(t)^2);已知l=10；R2=5;R1=10;a=1；积分下限为0；求上限值.

Answer 1

令f(t)=√(x^2+y^2+z^2)*x,易知f(t)是以2*π为周期的周期函数
再令变上限积分函数：F(x)=∫{0,x}f(t)dt,x≥0
可以证明F(x)是以2*π为周期的周期函数(见附录),F(x)图像如下图所示

 
因此,若F(b)=∫{0,b}f(t)dt=10且该方程有解,则在b≥0时有无穷多解
而对一个具体问题,若积分上限有范围,求出该范围内的解即可
若没有范围限制,求出一个周期内的所有解,再加上周期的整数倍即可
从上图可见,b在[0,1]与[2,3]均有一根,即一个周期内有两个实根
 
函数m文件(存为ni.m)：
----------------------------------------------------------------------------------------
function I=ni(b)
%用辛普森公式求数值积分∫{0,b}f(t)dt
for i=1:length(b)
    I(i)=quad(@integrand,0,b(i));
end
 
function f=integrand(t)
%建立被积函数f(t)
a=1;
R1=10;
R2=5;
x=R2*cos(t)/a+sqrt(a^2+1)*cos(t).*sin(t)./(a*(R1^2-R2^2*cos(t).^2));
y=R2*cos(t);
z=-R2*sin(t);
f=sqrt(x.^2+y.^2+z.^2).*x;
----------------------------------------------------------------------------------------
输入：
clear
fz=@(b)(ni(b)-10);   %建立方程I=10
b1=fzero(fz,[0,1])     %有根区间取[0,1],以下类似
b2=fzero(fz,[2,3])
b3=fzero(fz,[6,7])
b4=fzero(fz,[9,10])
 
输出：
b1=0.288539440045146
b2=2.853053213392633
b3=6.571724770845351
b4=9.136238520186859
----------------------------------------------------------------------------------------
可见b3≈b1+2*π,b2≈b4+2*π
 
附录　证明F(x)是以2*π为周期的周期函数
∵∫{0,2*π}f(t)dt=∫{0,π}f(t)dt+∫{π,2*π}f(t)dt
                       =∫{0,π}[f(t)+f(-t)]dt                                    对第2个积分换元并整理得
                       =∫{0,π/2}[f(t)+f(-t)]dt+∫{π/2,π}[f(t)+f(-t)]dt
                       =∫{0,π/2}[f(t)+f(π-t)+f(-t)+f(t-π)]dt            对第2个积分换元并整理得
                       =0                                                           由于f(t)+f(π-t)=0,f(-t)+f(t-π)=0
 
∴F(x+2*π)=∫{0,x+2*π}f(t)dt
                 =∫{0,x }f(t)dt+∫{x,x+2*π}f(t)dt
                 =∫{0,x }f(t)dt+∫{0,2*π}f(t)dt                  由于f(t)是以2*π为周期的周期函数
                 =∫{0,x }f(t)dt+0
                 =F(x)
即F(x)是以2*π为周期的周期函数