声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3526|回复: 6

[编程技巧] ode45函数调用问题请教。

[复制链接]
发表于 2009-9-4 23:46 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
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                  
请高手指点,该怎么解决?
回复
分享到:

使用道具 举报

发表于 2009-9-5 00:11 | 显示全部楼层
没空细看, 直觉LZ的函数有问题!?
函数里头怎有符号(symbolic)?
先看看help ode45!
 楼主| 发表于 2009-9-5 00:24 | 显示全部楼层


程序从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)
可以顺利执行,但感觉这种办法效率很低,纯手动的,有没有好点的解决办法呢?
谢谢!
 楼主| 发表于 2009-9-13 09:38 | 显示全部楼层
最近又有个新的问题:微分方程必须要先写成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);
发表于 2009-9-13 11:00 | 显示全部楼层
函数给定的方法有很多种, 如m文件/Inline/arrayfun/匿名函数(Anonymous Function)/内嵌函数(Nested Function), 愈新版本方式愈多!
其中许多个人也不熟, ODE/PDE这一块来此前没用过, LZ可搜下, 好像有人提过比较!
 楼主| 发表于 2009-9-13 11:05 | 显示全部楼层
原帖由 ChaChing 于 2009-9-13 11:00 发表
函数给定的方法有很多种, 如m文件/Inline/arrayfun/匿名函数(Anonymous Function)/内嵌函数(Nested Function), 愈新版本方式愈多!
其中许多个人也不熟, ODE/PDE这一块来此前没用过, LZ可搜下, 好像有人提过比较!


谢谢你,之前搜过关于ODE的帖子,但一直也没有发现有什么帮助。呵呵
 楼主| 发表于 2009-9-15 16:09 | 显示全部楼层
补充个问题:
形如dy/dx=f(x),这样的微分方程,怎么用ode45来求解?是否一定需要个y=g(x)的方程,才能用ode45求解?谢谢!
ode45求解的都是形如dy/dx=f(x,y),这样的方程?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-6 00:38 , Processed in 0.069749 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表