ode45函数调用问题请教。

quakefans

function程序如下：
function [y,u]=du(y,u)
syms y s t dta dtac reo pd;
s1='s^3+t*s^2/(dta^(2/3))-1=0';
t=5;
dta=200;
reo=458.29;
subs(s1);                %
s1=ans;
v=solve(s1,'s');
for k=1:length(v)
    idx(k) = isreal(v(k,1));
end
z=v(idx);
z
pd='((3/4)*reo)^(1/3)-t';
subs(pd);                %
pd=ans;
pd
if pd>0
   dtac='(1/(s^3/2))*(((3/4)*reo)^(1/3)-t)^(3/2)';
   subs(dtac,s,z);
   dtac=ans;
   subs(dtac);
   dtac=ans;
   dtac
else
   dtac=21;
   dtac
end
eddv='-0.5+0.5*(1+0.64*y^2*(1-s^3*y/dta)^2*(1-exp((-y)/26*(1-s^3*y/dta)^(1/2)*(1-dtac/dta)))^2)^(1/2)';
dta=200;
subs(eddv);
eddv=ans;
subs(eddv,s,z);
eddv=ans;
eddv
du=(1-s^3*y/dta)/(1+eddv);
dta=200;
subs(du);
du=ans;
subs(du,s,z);
du=ans;
du

在命令窗口中输入[y,u]=ode45('du',[0 1],1)
后出现如下错误：
??? Error using ==> odearguments
Inputs to odearguments must be floats, namely single or double.
Error in ==> odearguments at 136
  dataType = superiorfloat(t0,y0,f0);
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

错误信息说变量数据类型不对，我查了下数据类型，如下：
Name      Size            Bytes  Class      Attributes
  ans       1x1               116  sym                  
  dta       1x1                 8  double               
  dtac      1x1                 8  double               
  du        1x1               116  sym                  
  eddv      1x1               116  sym                  
  idx       1x3                 3  logical              
  k         1x1                 8  double               
  pd        1x1                 8  double               
  reo       1x1                 8  double               
  s         1x1                58  sym                  
  s1        1x1               116  sym                  
  t         1x1                 8  double               
  v         3x1               236  sym                  
  y         1x1                58  sym                  
  z         1x1               116  sym                  
请高手指点，该怎么解决？

ChaChing

没空细看, 直觉LZ的函数有问题!?
函数里头怎有符号(symbolic)?
先看看help ode45!

quakefans

原帖由 ChaChing 于 2009-9-5 00:11 发表

                               
登录/注册后可看大图

没空细看, 直觉LZ的函数有问题!?
函数里头怎有符号(symbolic)?
先看看help ode45!

程序从function [y,u]=du(y,u)行之后，就是经过一些运算求出函数表达式du，运算过程中，有符号代换的运算。

刚才试下了，比较笨的办法，就是把运算出的du的表达式，直接写成函数：
function [y,u]=du(y,u)
du=-((y*(200^(2/3)/(14400*((2159^(1/2)*8640^(1/2))/8640 + 4319/8640)^(1/3)) - 200^(1/3)/120 + ((2159^(1/2)*8640^(1/2))/8640 + 4319/8640)^(1/3))^3)/200 - 1)/(0.5*(0.64*y^2*(1/exp((217807536381482829*y*(1 - (y*(200^(2/3)/(14400*((2159^(1/2)*8640^(1/2))/8640 + 4319/8640)^(1/3)) - 200^(1/3)/120 + ((2159^(1/2)*8640^(1/2))/8640 + 4319/8640)^(1/3))^3)/200)^(1/2))/5854679515581644800) - 1)^2*((y*(200^(2/3)/(14400*((2159^(1/2)*8640^(1/2))/8640 + 4319/8640)^(1/3)) - 200^(1/3)/120 + ((2159^(1/2)*8640^(1/2))/8640 + 4319/8640)^(1/3))^3)/200 - 1)^2 + 1)^(1/2) + 0.5);

然后再执行 [y,u]=ode45('du',[0:0.01:1],1)
可以顺利执行，但感觉这种办法效率很低，纯手动的，有没有好点的解决办法呢？
谢谢！

quakefans

最近又有个新的问题：微分方程必须要先写成function然后再在指令窗口用ode45求解吗？直接在指令窗口写函数，给出函数表达式后，紧接着用ode45命令貌似也可以求解。这里有点不太明白，请指教，谢谢！
syms y s t dta dtac reo pd;
s1='s^3+t*s^2/(dta^(2/3))-1=0';
t=0;
dta=10;
reo=458.29;
subs(s1);                %
s1=ans;
v=solve(s1,'s');
for k=1:length(v)
    idx(k) = isreal(v(k,1));
end
z=v(idx);
%z=vpa(z,32);
z
pd='((3/4)*reo)^(1/3)-t';
subs(pd);                %
pd=ans;
pd
if pd>0
   dtac='(1/(s^3/2))*(((3/4)*reo)^(1/3)-t)^(3/2)';
   subs(dtac,s,z);
   dtac=ans;
   subs(dtac);
   dtac=ans;
   dtac
else
   dtac=21;
   dtac
end
eddv='-0.5+0.5*(1+0.64*y^2*(1-s^3*y/dta)^2*(1-exp((-y)/26*(1-s^3*y/dta)^(1/2)*(1-dtac/dta)))^2)^(1/2)';
%dta=200;
subs(eddv);
eddv=ans;
subs(eddv,s,z);
eddv=ans;
%eddv=vpa(eddv,8);
eddv
du=(1-s^3*y/dta)/(1+eddv);
%dta=200;
%subs(du);
%du=ans;
subs(du,s,z);
du=ans;
du                                         %到这里，函数的表达式du已经计算出来，下面紧接着就用ode45命令。
[y,u]=ode45('du',[0 dta],1);

ChaChing

函数给定的方法有很多种, 如m文件/Inline/arrayfun/匿名函数(Anonymous Function)/内嵌函数(Nested Function), 愈新版本方式愈多!
其中许多个人也不熟, ODE/PDE这一块来此前没用过, LZ可搜下, 好像有人提过比较!

quakefans

原帖由 ChaChing 于 2009-9-13 11:00 发表

                               
登录/注册后可看大图

函数给定的方法有很多种, 如m文件/Inline/arrayfun/匿名函数(Anonymous Function)/内嵌函数(Nested Function), 愈新版本方式愈多!
其中许多个人也不熟, ODE/PDE这一块来此前没用过, LZ可搜下, 好像有人提过比较!

谢谢你，之前搜过关于ODE的帖子，但一直也没有发现有什么帮助。呵呵

quakefans

补充个问题：
形如dy/dx=f(x)，这样的微分方程，怎么用ode45来求解？是否一定需要个y=g(x)的方程，才能用ode45求解？谢谢！
ode45求解的都是形如dy/dx=f(x,y)，这样的方程？