西西想幸福 发表于 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 # 西西想幸福 的帖子

是啊,这个并不算难
页: 1 [2] 3
查看完整版本: 帮个忙吧,用ode45解微分方程组,刚开始学,请指点一二!!