帮个忙吧，用ode45解微分方程组，刚开始学，请指点一二！！

西西想幸福

x(3)不是x的二阶导数啊！！我也是看了好多ode的例子然后编出的这个方程！！那要不我试试用矩阵的形式再编一下？那样会有不一样么？

西西想幸福

我又把方程改成这样的了！！可是还是不行啊！！这回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

西西想幸福

问一下，要是想让计算停止要用什么语句？

321forever

倒是看懂了lz的方程，有点像是控制里面的空间状态方程啊，状态量是x1,x2,x3,tg,ts.
感觉编程也没有什么问题，我换了ode15s试了下，倒是出图了，lz可以看下图，但我觉得图很奇怪，lz再看看能不能再简化下方程
[t,x]=ode15s(@modfun,tspan,x0);  

西西想幸福

是不是直接把那个ode45换成这个就行啊？是像下面这个方程这样么？
clc
clear
tspan=[0:0.01:1000];                                                                          %定义变量求解区间
x0=[1.4*10^-6,1.5*10^-6,0,650,298];                                                            %设置初值，并设在z=0处x1的值是1.4*10^-6，x2的值是1.5*10^-6,tg是650k
options=odeset('RelTol',1e-3);                                                           %设置相对误差
[t,x]=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
>>

那里错了？

西西想幸福

clc
clear
tspan=[0:0.01:1000];                                                                        
x0=[1.4*10^-6,1.5*10^-6,0,650,298];                                                            
options=odeset('RelTol',1e-3);                                                        
[t,x]=ode15s(@modfun,tspan,x0);                                                  
plot(t,x(:,1),'k:');
hold on;
plot(t,x(:,3),'k:');
set(gca,'Fontsize',12);
xlabel('\itt','Fontsize',16);


有些是我后改的！！原来的程序只能出来一条虚线！！应该是不对！！

西西想幸福

要不我明天再从新化简一下方程！！再试试！！那后面这个错误提示你的有么？这个是什么意思？

321forever

回复 22 # 西西想幸福 的帖子

我也出现这个warning，我觉得方程很复杂，简化下再看看

西西想幸福

恩，好的！！嘻嘻！！今天就化简下！！谢谢！！

西西想幸福

今天又从新化简了一下方程，更正了一些错误，可是结果还是这样啊！！怎么办啊！！这个warning到底是什么意思呢？是我的方程有问题么？

segev

楼主应该把数学方程式用公式编辑器列出来，直接看代码看不明白

西西想幸福

这是我的方程截图

是关于冶金动力学的一个方程！！大家再帮忙看下！！有没有人是学冶金的啊？因为我觉得我的方程参数可能代错了！！有没有人懂啊？帮帮忙吧！！

segev

回复 27 # 西西想幸福 的帖子

并非此专业，看不懂这个方程，但是你描述的问题貌似可能是遇到了一个刚性问题，用ode23t试试，或者自己写一个算法也可以

西西想幸福

自己写算法是指自己编一个龙格-库塔的程序么？

segev

回复 29 # 西西想幸福 的帖子

是啊，这个并不算难