MATLAB 求积分时无法输出数值解,一直是积分表达式syms t x0=1753000;y0=0;z0=0;vx0=0;vy0=1700;vz0=0;x1=1738000*cosd(10);y1=1738000*sind(10);z1=0;vx1=0;vy1=0;vz1=0;u=4.9*10^12 ;r0=1753000;a0=1.6243;b0=0;c0=0;a1=(-6*(a0*t^2+(vx1+3*vx0)*t-

MATLAB 求积分时无法输出数值解,一直是积分表达式syms t x0=1753000;y0=0;z0=0;vx0=0;vy0=1700;vz0=0;x1=1738000*cosd(10);y1=1738000*sind(10);z1=0;vx1=0;vy1=0;vz1=0;u=4.9*10^12 ;r0=1753000;a0=1.6243;b0=0;c0=0;a1=(-6*(a0*t^2+(vx1+3*vx0)*t-
MATLAB 求积分时无法输出数值解,一直是积分表达式
syms t
x0=1753000;
y0=0;
z0=0;
vx0=0;
vy0=1700;
vz0=0;
x1=1738000*cosd(10);
y1=1738000*sind(10);
z1=0;
vx1=0;
vy1=0;
vz1=0;
u=4.9*10^12 ;
r0=1753000;
a0=1.6243;
b0=0;
c0=0;
a1=(-6*(a0*t^2+(vx1+3*vx0)*t-4*(x1-x0)))/t^3;
b1=(-6*(b0*t^2+(vy1+3*vy0)*t-4*(x1-x0)))/t^3;
c1=(-6*(c0*t^2+(vz1+3*vz0)*t-4*(x1-x0)))/t^3;
a2=(6*(a0*t^2+2*(vx1+2*vx0)*t-6*(x1-x0)))/t^4;
b2=(6*(b0*t^2+2*(vy1+2*vy0)*t-6*(y1-y0)))/t^4;
c2=(6*(c0*t^2+2*(vz1+2*vz0)*t-6*(z1-z0)))/t^4;
ax=a0+a1*t+a2*t^2;
ay=b0+b1*t+b2*t^2;
az=c0+c1*t+c2*t^2;
x=x0+vx0*t+a0*t^2/2+a1*t^3/6+a2*t^4/12;
y=y0+vy0*t+b0*t^2/2+b1*t^3/6+b2*t^4/12;
z=z0+vz0*t+c0*t^2/2+c1*t^3/6+c2*t^4/12;
r=sqrt(x^2+y^2+z^2)
gx=(-u*x)/r^3;
gy=(-u*y)/r^3;
gz=(-u*z)/r^3;
apx=ax-gx;
apy=ay-gy;
apz=az-gz;
ap=sqrt(apx^2+apy^2+apz^2);
V=int(ap,t,0,200);
试过eval函数和vpa函数,但是V一直显示为int(……,t,0,200)而不能得到具体数值,

MATLAB 求积分时无法输出数值解,一直是积分表达式syms t x0=1753000;y0=0;z0=0;vx0=0;vy0=1700;vz0=0;x1=1738000*cosd(10);y1=1738000*sind(10);z1=0;vx1=0;vy1=0;vz1=0;u=4.9*10^12 ;r0=1753000;a0=1.6243;b0=0;c0=0;a1=(-6*(a0*t^2+(vx1+3*vx0)*t-
syms t 
x0=1753000;
y0=0;
z0=0;
vx0=0;
vy0=1700;
vz0=0;
x1=1738000*cosd(10);
y1=1738000*sind(10);
z1=0;
vx1=0;
vy1=0;
vz1=0;
u=4.9*10^12 ;
r0=1753000;
a0=1.6243;
b0=0;
c0=0;
a1=(-6*(a0*t^2+(vx1+3*vx0)*t-4*(x1-x0)))/t^3;
b1=(-6*(b0*t^2+(vy1+3*vy0)*t-4*(x1-x0)))/t^3;
c1=(-6*(c0*t^2+(vz1+3*vz0)*t-4*(x1-x0)))/t^3;
a2=(6*(a0*t^2+2*(vx1+2*vx0)*t-6*(x1-x0)))/t^4;
b2=(6*(b0*t^2+2*(vy1+2*vy0)*t-6*(y1-y0)))/t^4;
c2=(6*(c0*t^2+2*(vz1+2*vz0)*t-6*(z1-z0)))/t^4;
ax=a0+a1*t+a2*t^2;
ay=b0+b1*t+b2*t^2;
az=c0+c1*t+c2*t^2;
x=x0+vx0*t+a0*t^2/2+a1*t^3/6+a2*t^4/12;
y=y0+vy0*t+b0*t^2/2+b1*t^3/6+b2*t^4/12;
z=z0+vz0*t+c0*t^2/2+c1*t^3/6+c2*t^4/12;
r=sqrt(x^2+y^2+z^2);
gx=(-u*x)/r^3;
gy=(-u*y)/r^3;
gz=(-u*z)/r^3;
apx=ax-gx;
apy=ay-gy;
apz=az-gz;
ap=sqrt(apx^2+apy^2+apz^2);

fun=inline(vectorize(char(eval(ap))),'t');
tt=0:200;
yy=fun(tt);
plot(tt,yy);
V=quad(fun,1,200);


由于函数太复杂,不能得到精确的解析积分
这里可以完全用数值办法
fun=inline(vectorize(char(eval(ap))),'t');

将原来ap的内容转为inline函数fun(t)


我们可以画出其在0~200范围内的图像
当tt趋向于0时,yy趋向于无穷大,所以函数的0点奇点


[0 200]区间内的积分是发散的
所以这里计算了[1 200]区间内的积分
V=quad(fun,1,200);