ode45函数调用问题请教。
function程序如下：function =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
在命令窗口中输入=ode45('du',,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 BytesClass Attributes
ans 1x1 116sym
dta 1x1 8double
dtac 1x1 8double
du 1x1 116sym
eddv 1x1 116sym
idx 1x3 3logical
k 1x1 8double
pd 1x1 8double
reo 1x1 8double
s 1x1 58sym
s1 1x1 116sym
t 1x1 8double
v 3x1 236sym
y 1x1 58sym
z 1x1 116sym
请高手指点，该怎么解决？ 没空细看, 直觉LZ的函数有问题!?
函数里头怎有符号(symbolic)?
先看看help ode45! 原帖由 ChaChing 于 2009-9-5 00:11 发表 http://www.chinavib.com/forum/images/common/back.gif
没空细看, 直觉LZ的函数有问题!?
函数里头怎有符号(symbolic)?
先看看help ode45!
程序从function =du(y,u)行之后，就是经过一些运算求出函数表达式du，运算过程中，有符号代换的运算。
刚才试下了，比较笨的办法，就是把运算出的du的表达式，直接写成函数：
function =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);
然后再执行 =ode45('du',,1)
可以顺利执行，但感觉这种办法效率很低，纯手动的，有没有好点的解决办法呢？
谢谢！ 最近又有个新的问题：微分方程必须要先写成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命令。
=ode45('du',,1); 函数给定的方法有很多种, 如m文件/Inline/arrayfun/匿名函数(Anonymous Function)/内嵌函数(Nested Function), 愈新版本方式愈多!
其中许多个人也不熟, ODE/PDE这一块来此前没用过, LZ可搜下, 好像有人提过比较! 原帖由 ChaChing 于 2009-9-13 11:00 发表 http://www.chinavib.com/forum/images/common/back.gif
函数给定的方法有很多种, 如m文件/Inline/arrayfun/匿名函数(Anonymous Function)/内嵌函数(Nested Function), 愈新版本方式愈多!
其中许多个人也不熟, ODE/PDE这一块来此前没用过, LZ可搜下, 好像有人提过比较!
谢谢你，之前搜过关于ODE的帖子，但一直也没有发现有什么帮助。呵呵 补充个问题：
形如dy/dx=f(x)，这样的微分方程，怎么用ode45来求解？是否一定需要个y=g(x)的方程，才能用ode45求解？谢谢！
ode45求解的都是形如dy/dx=f(x,y)，这样的方程？
页:
[1]