ode求解的m函数定义问题

sindy209

function dy=mylu(t,Y)
a=Y(1);
b=Y(2);
c=Y(3);
dY(1)=-1.5-b*(5*a^4+2*c+1);
dY(2)=-a^5-2*a*c-a;
dY(3)=-2*a*b;
我这样定义有错？
为什么调用的时候会出现Input argument "Y" is undefined.
公式如 Doc11.doc (17.5 KB, 下载次数: 14)

2008-4-12 20:25 上传

点击文件名下载附件

高手帮忙看下，谢谢！

[ 本帖最后由 无水1324 于 2008-4-14 19:04 编辑 ]

咕噜噜

你本来就没定义，不过你最好把Y（1）等直接带进方程中，不要设置成a等符号
不过你单独运行子程序肯定会告诉你没定义，这是正常的

sindy209

子程序：
function dy=mylu(t,Y)
dY=zeros(3,3);
x=zeros(3,3);

Y=-1.5*x(1)-x(4)*(5*x(1)^4+2*x(5)+1);
dY(1)=-1.5-x(4)*(5*x(1)^4+2*x(5)+1);
dY(2)=-x(1)^5-2*x(1)*x(5)-x(1);
dY(3)=-2*x(1)*x(4);
q1=x(1)^2*(x(1)^4+2*x(5)-1)-0.001*x(4);
q2=x(1)^2*(x(1)^4+2*x(5)-1)-0.001*x(4)+(x(2)+1.5*x(1)+x(4)*x(1)*(x(1)^4+2*x(5)+1))^2*((1+dY(1)^2)*(2*x(5)+1)+(x(1)^2+x(2)^2)+(dY(1)*X(1)^2)^2);
主程序：
怎么使用来调用子程序中的dy(1),dy(2),dy(3),q1,q2呢？
谢谢！
主程序：
clear all

lubang=inline('[x(2)+x(1)^2*sin(2*x(1))+2*x(3)*sin(x(1))+0.5*x(1)*sin(t);gama+5*x(1)+(6*x(2)^2)/(1+2*x(1)^2)+(1-cos(x(2)))*x(3)+0.5*x(1)*cos(t);-x(3)+x(1)^2+0.5;0.5*q2;-x(5)+x(1)^4+0.625]','t','x');

gama=-2*(xa(:,2)-Y)^2+dy(1)*xa(:,2)+0.5*dY(2)*q+dY(3)*xa(:,5)-xa(:,4).*(xa(:,2)-Y)^2*(1+dY(1)^2)*(2*xa(:,5)+1)-xa(:,4).*(xa(:,2)-Y)^2*(((xa(:,1).^2+xa(:,2).^2).^2)^2+(dY(1)*xa(:,1).^2)^2);
options=odeset('RelTol',1e-4,'AbsTol',1e-4);
[t1,xa]=ode45(lb,[0 2],[1,1,1,1,1,1],options);
figure(1)
plot(t2,xa(:,1),'k:')

[ 本帖最后由 咕噜噜 于 2008-4-14 15:55 编辑 ]

咕噜噜

好乱:@L :@L 看看
子程序和主程序都没看明白

[ 本帖最后由 咕噜噜 于 2008-4-14 15:57 编辑 ]

sindy209

子程序1：
function dy=mylu1(t,x)
dy=zeros(3,1);
dy(1)=-1.5-x(2)*(5*x(1).^4+2*x(3)+1);
dy(2)=-x(1).^5-2*x(1)*x(3)-x(1);
dy(3)=-2*x(1)*x(2);

主程序：
clear all
x=zeros(5,5);
[t,dy]=ode45('mylu1',[0 2],[-9.5,-4,-2]);
Y=-1.5*x(1)-x(2)*(5*x(1).^4+2*x(3)+1);
q=x(1).^2*(x(1).^4+2*x(3)-1)-0.001*x(2)+(x(4)+1.5*x(1)+x(2)*x(1)*(x(1)^4+2*x(3)+1)).^2*((1+dy(1)^2)*(2*x(3)+1)+(x(1).^2+x(4).^2)+(dy(1)*x(1).^2).^2);


gama=-2*(x(4)-Y).^2+dy(1)*x(4)+0.5*dy(2)*q+dy(3)*x(3)-x(2)*(x(4)-Y)*(1+dy(1)^2)*(2*x(3)+1)-x(4)*(x(4)-Y)*(((x(1).^2+x(4).^2).^2)^2+(dy(1)*x(1).^2)^2);

lubang=inline('[x(4)+x(1).^2*sin(2*x(1))+2*x(5)*sin(x(1))+0.5*x(1)*sin(t);0.5*q;-x(3)+x(1).^4+0.625;gama+5*x(1)+(6*x(4).^2)/(1+2*x(1).^2)+(1-cos(x(4)))*x(5)+0.5*x(1)*cos(t);-x(5)+x(1).^2+0.5]','t','x');
options=odeset('RelTol',1e-4,'AbsTol',1e-4);
[t1,x]=ode45(lubang,[0 2],[1,1,1,1,1,1],options);

figure(1)
plot(t2,x(:,1),'k:')

运行出现：
Warning: Failure at t=5.935321e-003.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.387779e-017) at time t.
> In ode45 at 378
Warning: Failure at t=5.935321e-003.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.387779e-017) at time t.
> In ode45 at 378
??? Undefined command/function 'x'.
不是定义x=zeros(5,5)了？
步长应该怎么减呢？


帮忙看下，弄了好久都……

在这里特别感谢octopussheng 和无水，咕噜噜的帮忙！谢谢大家了！
在此谢谢啦！

[ 本帖最后由 sindy209 于 2008-4-14 22:57 编辑 ]

无水1324

function dy=mylu1(t,x)
dy=[-1.5-x(2)*(5*x(1).^4+2*x(3)+1);
    -x(1).^5-2*x(1)*x(3)-x(1);
    -2*x(1)*x(2)];

但是你的程序有问题求解之后解瞬间变得很大

sindy209

楼主能说的更明白点嘛！
是mylu1.m这个程序有错？
谢谢！

无水1324

对
按照我6楼的，使可以运行计算出来，但是解可能不是你想要的

咕噜噜

他的公式传上来，以便大家参考

咕噜噜

所以6楼的程序是正确的
但是这么求解会不会有问题，我是说公式

[ 本帖最后由 无水1324 于 2008-4-14 19:10 编辑 ]

octopussheng

最初的问题很简单
dY(1)=-1.5-b*(5*a^4+2*c+1);
dY(2)=-a^5-2*a*c-a;
dY(3)=-2*a*b;
把Y改成y就可以了！

sindy209

修改后的程序如5楼
可是还是有错:@Q :'( ，麻烦大家在帮我看下！
谢谢！

无水1324

同学你看一下六楼的程序怎么样？

sindy209

按照6楼的方法做，会出现这样的错误
Warning: Failure at t=5.935321e-003.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.387779e-017) at time t.
> In ode45 at 378
??? Error using ==> inlineeval
Error in inline expression ==> [x(4)+x(1).^2*sin(2*x(1))+2*x(5)*sin(x(1))+0.5*x(1)*sin(t);0.5*q;-x(3)+x(1).^4+0.625;gama+5*x(1)+(6*x(4).^2)/(1+2*x(1).^2)+(1-cos(x(4)))*x(5)+0.5*x(1)*cos(t);-x(5)+x(1).^2+0.5]
??? Error using ==> eval
Undefined function or variable 'q'.

Error in ==> inline.feval at 34
        INLINE_OUT_ = inlineeval(INLINE_INPUTS_, INLINE_OBJ_.inputExpr, INLINE_OBJ_.expr);

Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Y也是表达式，为什么就没有错误，而q也是表达式，为什么就错了呢？

咕噜噜

这种问题没有遇到过，应该是没问题啊，是不是那里有地方不一致
你不写你新的程序看不到拉
显示q没有定义