五散人 发表于 2014-4-30 14:34

简单微分方程组求解


大家好,我对某个简单阀列些了微分方程组,使用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;endendpump 函数 function y= pump(x)%UNTITLED5 Summary ofthis function goes here%   Detailed explanation goes hereif x>=1.8e7   y=0;else   y=6.67e-3;end
end

lihaitao123 发表于 2014-5-8 10:13

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

五散人 发表于 2014-5-8 20:58

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

gghhjj 发表于 2014-5-18 01:35

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

Kevin_HIT 发表于 2014-10-9 21:13

学习了学习了

Kevin_HIT 发表于 2014-10-13 19:42

学习了学习了

孟铭恩 发表于 2014-10-20 21:42

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

zhngyongshun 发表于 2015-1-16 10:40

问题解决了么
页: [1]
查看完整版本: 简单微分方程组求解