西西想幸福
发表于 2012-3-20 17:16
x(3)不是x的二阶导数啊!!我也是看了好多ode的例子然后编出的这个方程!!那要不我试试用矩阵的形式再编一下?那样会有不一样么?
西西想幸福
发表于 2012-3-20 17:34
我又把方程改成这样的了!!可是还是不行啊!!这回conmand窗口里什么也没有了!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dx=modfun(t,x)
rc=(0.125-3.78*x(3))^1/3;
d1=1.467*10^-6*x(4)^1.75;
d2=3.828*10^-7*x(4)^1.75;
k1=0.225*exp(-179.14/x(4));
k2=0.650*exp(-342.4324/x(4));
r1=-(4*pi*rc^2*(0.529-x(1)))/(1/k1+rc/d1-4*rc^2/d1);
r2=-(4*pi*rc^2*(0.347-x(2)))/(1/k2+rc/d2-4*rc^2/d2);
cpg=104.1273+0.9806*x(4)-363688*x(4)^-2-314.5*x(1)-655.2*x(2)-0.3224*x(1)*x(4)-2.8364*x(4)*x(2)+24000*x(4)^-2*x(1)+9372000*x(4)^-2*x(2);
gms=2533.3*(-4.4864*rc^3+0.0069)/(-0.0342*rc^3+0.0069);
cps=(3.7125+0.006*x(5)*rc^3+79.9140*rc^3)/(0.3875-0.3*rc^3);
h1=999.8621-3.79*x(5)+0.0011*x(5)^2+1.2*10^4*x(5)^-1;
h2=840.3932-0.79*x(5)+5.333*10^-4*x(5)^2-1.95*10^5*x(5)^-1;
dx=[-0.99/1.98*10^8*r1;...
-0.99/1.98*10^8*r2;...
-0.99/4.56*10^4*(r1+r2);...
14.07*(x(5)-x(4))/cpg;...
0.99/gms/cps*(pi*4*10^4*(x(5)-x(4))-h1*r1-h2*r2)];
dx=dx(:);
return
西西想幸福
发表于 2012-3-20 17:34
问一下,要是想让计算停止要用什么语句?
321forever
发表于 2012-3-20 19:45
倒是看懂了lz的方程,有点像是控制里面的空间状态方程啊,状态量是x1,x2,x3,tg,ts.
感觉编程也没有什么问题,我换了ode15s试了下,倒是出图了,lz可以看下图,但我觉得图很奇怪,lz再看看能不能再简化下方程
=ode15s(@modfun,tspan,x0);
西西想幸福
发表于 2012-3-20 21:10
是不是直接把那个ode45换成这个就行啊?是像下面这个方程这样么?
clc
clear
tspan=; %定义变量求解区间
x0=; %设置初值,并设在z=0处x1的值是1.4*10^-6,x2的值是1.5*10^-6,tg是650k
options=odeset('RelTol',1e-3); %设置相对误差
=ode15s(@modfun,tspan,x0); %调用ode45函数解微分方程
figure;
plot(t,x(:,1),'k:');
hold on;
plot(t,x(:,3),'k:');
set(gca,'Fontsize',12);
xlabel('\itt','Fontsize',16);
可是为什么会出现这个?
Warning: Failure at t=1.684686e+000.Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (3.552714e-015) at time t.
> In ode15s at 753
In modfun1 at 6
>>
那里错了?
西西想幸福
发表于 2012-3-20 21:14
clc
clear
tspan=0.01:1000];
x0=;
options=odeset('RelTol',1e-3);
=ode15s(@modfun,tspan,x0);
plot(t,x(:,1),'k:');
hold on;
plot(t,x(:,3),'k:');
set(gca,'Fontsize',12);
xlabel('\itt','Fontsize',16);
有些是我后改的!!原来的程序只能出来一条虚线!!应该是不对!!
西西想幸福
发表于 2012-3-20 21:26
要不我明天再从新化简一下方程!!再试试!!那后面这个错误提示你的有么?这个是什么意思?
321forever
发表于 2012-3-20 21:42
回复 22 # 西西想幸福 的帖子
我也出现这个warning,我觉得方程很复杂,简化下再看看
西西想幸福
发表于 2012-3-21 08:17
恩,好的!!嘻嘻!!今天就化简下!!谢谢!!
西西想幸福
发表于 2012-3-22 17:00
今天又从新化简了一下方程,更正了一些错误,可是结果还是这样啊!!怎么办啊!!这个warning到底是什么意思呢?是我的方程有问题么?
segev
发表于 2012-4-3 15:27
楼主应该把数学方程式用公式编辑器列出来,直接看代码看不明白
西西想幸福
发表于 2012-4-5 09:13
本帖最后由 西西想幸福 于 2012-4-5 09:14 编辑
是关于冶金动力学的一个方程!!大家再帮忙看下!!有没有人是学冶金的啊?因为我觉得我的方程参数可能代错了!!有没有人懂啊?帮帮忙吧!!
segev
发表于 2012-4-5 16:15
回复 27 # 西西想幸福 的帖子
并非此专业,看不懂这个方程,但是你描述的问题貌似可能是遇到了一个刚性问题,用ode23t试试,或者自己写一个算法也可以
西西想幸福
发表于 2012-4-6 08:35
自己写算法是指自己编一个龙格-库塔的程序么?
segev
发表于 2012-4-6 14:57
回复 29 # 西西想幸福 的帖子
是啊,这个并不算难