简单微分方程组求解

五散人

大家好，我对某个简单阀列些了微分方程组，使用ODE45、ode23s、Ode15s等解法，均失败。ode15s提示 Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-019) at time t. 而ode23s，一直处于busy。下面贴出公式，和matlab代码，请高手不吝赐教。

参数表：

pi=3.14;

alpha=pi/6;%半锥角

rou=1000;%密度

Beta=7e8;%体积弹性模量

V01=8.792e-3;%管道容积

V02=1.4e-5;%阀芯后部容积

m=0.2;%阀芯质量

C=0.8;%锥阀流量系数

d=2.6e-2;%阀芯平均直径

K1=2e4;%弹簧刚度

x0=3e-3;%弹簧预压缩量

Ks=C*pi*d*sin(2*alpha);%液动力系数

Cd=C*pi*d*sin(alpha);%锥阀流量计算

A2=6.15e-4;%阀芯面积

Q0=6.67e-3;%泵流量

B=0.02;%粘度系数

f=30;%摩擦力

d1=1.5e-3;%阻尼孔1直径

l1=1.2e-2; %阻尼孔1长度

d2=2.5e-3;%阻尼孔2直径

l2=8e-3; %阻尼孔2长度

Ea1=1e8;%壁面刚度

Ea2=0.8e7;%阀座刚度

G1=0.1351*power(d1^8/l1,1/3);%阻尼孔1流导

G2=0.1351*power(d2^8/l2,1/3);%阻尼孔2流导

dy=zeros(4,1);

dy(1)=y(2);

dy(2)=(y(3)*A2-y(4)*A2-K1*(y(1)+x0)-B*y(2)-Ks*y(1)*y(3)-f*sign(y(2))-Ea1*touch(y(1))-Ea2*touch(y(1)))/m;

dy(3)=(pump(y(3))-Cd*y(1)*sqrt(2*y(3)/rou)-G1*power((y(3)-y(4)),2/3)-A2*y(2))*Beta/V01;

dy(4)=(G1*power((y(3)-y(4)),2/3)-G2*power(y(4),2/3)+A2*y(2))*Beta/V02;

touch函数

function y=touch(x)

if (x-0.011)>=0

   y=x-0.011;

elseif x<=0

   y=x;

else

   y=0;

end

end

pump 函数

function y= pump(x)

%UNTITLED5 Summary ofthis function goes here

%   Detailed explanation goes here

if x>=1.8e7

   y=0;

else

   y=6.67e-3;

end

end

lihaitao123

1不知道楼主的初值是否可以改变。2可以用newmark-beta或者可以调自定义步长例如RK-4 方法试一下。我之前碰到是不是楼主这种有实际意义的非线性方程，通过改变初值可以运行。

五散人

首先表示感谢。我尝试过改变初值，但因为是实际问题，受限于初始时位置、速度为0，两个压力相等的条件，没办法做太多调整。以前也用过列些微分方程组，求解系统动态。没想到这次中招了。

gghhjj

自己编写个算法程序，尝试跟中迭代过程中状态量的变化规律，看看是否能够找出问题所在

Kevin_HIT

学习了学习了

Kevin_HIT

学习了学习了

孟铭恩

Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-019) at time 在没有减小所允许最小步长的幅值，就不能满足整合公差？？   这啥意思   就看懂个加速度  都忘了

zhngyongshun

问题解决了么