你计算的公式应该是变量不只一个,然后两个数值差的比较大,导致经过相同的时间,两参数数值变化过大。比如说如下式子
Dx=@(t,x)[(x(1)-123*x(2)));
(x(2)-1234*x(1))];
[a,x]= ode45(Dx,,[0 36000],[80 3000]);
x=x(end,:)
t1=@(t)v(t)-150000; %%%%v(t)是计算t时刻储气室体积的函数
t_zong=fzero(t1,10000) %%%%求函数t1在t=10000附近的零点
这个瞎写的式子中,假设正确的结果是总时间t_zong=3600秒,即1个小时,x(1)代表温度,开始为80,1个小时后为160;x(2)代表一个储气室的体积,开始为2500立方米,1个小时后为15万立方米。那么输出结果应该为:
x=(160 150000)
t_zong=3600
这样的话,3600秒内x(1)才变化80,而x(2)会变化十几万,相差太大,程序很容易报上面你写的那个警告,无法运算下去。
那怎么解决呢?
很简单,把x(2)看成体积除以1000,算出来结果后,分别将x(2)中的各数值乘以1000就是你要算的数,即:
Dx=@(t,x)[(x(1)-123*x(2)*1000));
(x(2)*1000-1234*x(1))];
[a,x]= ode45(Dx,,[0 36000],[80 3]);
x=x(end,:)
t1=@(t)v(t)-150; %%%%v(t)是计算t时刻储气室体积的函数
t_zong=fzero(t1,10000) %%%%求函数t1在t=10000附近的零点
最后就应该能计算出结果了,结果应该是:
x=(160 150)
t_zong=3600
最后再把150乘上1000就是正确的体积了
具体除以多少是根据变化较小的那个确定的,在这个瞎写假设的情境中,温度变化是几十,那体积除以某个数后也应该是几十到几百。
把这个方法套到你的方程中,把某个数缩小或者放大,应该就能解决问题了。
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);
修改初始值