用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*π为周期的周期函数
2021-01-25 广告