用matlab 求定积分的上限b的值

令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*π为周期的周期函数