关于四阶五级Runge-Kutta变步长解微分方程的问题
这是我的程序%算变量E（t），在后面微分方程系数计算中使用
t=0:0.008:2; tn=t./(0.2+0.1555*0.8);
En=1.553174.*(tn./0.7).^1.9./(1+(tn./0.7).^1.9)./(1+(tn./1.173474).^21.9);
Ea=En.*3.0; E=Ea+0.06;
%采用ode45解微分方程（y=')
Vl=69.36; Qa=0; Pa=90; Pr=5;
if Pr>Pa
Dm=0; Da=1;
else
Dm=1; Da=0;
end% 循环是在这结束不?
%系数矩阵A和B
A11=-Da.*E; A12=0; A13=0; A14=Dm/0.8738;
A21=Da/0.001025.*E; A22=-(Da*0.01+0.0398)/0.001025; A23=-Da/0.001025; A24=0;
A31=0; A32=1/2.896; A33=-1/(0.8738.*2.896); A34=1/(0.8738.*2.896);
A41=Dm*E/4; A42=0; A43=-1/(0.8738*4); A44=0;
B1=-Dm.*(-E./5-0.06.*10); B2=-Da.*(-E./5-0.06.*10)/0.01;
B3=0; B4=-1./(4.*0.8738)-Dm./4.*(-E./5-0.06.*10+1);
%微分方程
Function dydt=equations(t,y)
A=[ A11 A12 A13 A14 ,
A21 A22 A23 A24 ,
A31 A32 A33 A34 ,
A41 A42 A43 A44 ];
B=;
dydt=A*y+B;
tspan=;
y0=;
options=odeset('Vl','Qa','Pa','Pr');
=ode45(@equations,tspan,y0,options);
plot(t,y(:,1))
为什么运行后说y是没有定义的？
ps：options=odeset('Vl','Qa','Pa','Pr')这个设置对吗？？谢谢！！
[ 本帖最后由 ChaChing 于 2009-4-17 08:41 编辑 ]
回复 #1 youurs 的帖子
??? Attempt to execute SCRIPT function as a function:D:\Program Files\MATLAB\R2007b\toolbox\matlab\lang\function.m
在红色加了end后
运行错误如上
Function dydt=equations(t,y)有错误 原帖由 youurs 于 2007-12-28 08:48 发表 http://www.chinavib.com/forum/images/common/back.gif
这是我的程序
%算变量E（t），在后面微分方程系数计算中使用
t=0:0.008:2;
tn=t./(0.2+0.1555*0.8);
En=1.553174.*(tn./0.7).^1.9./(1+(tn./0.7).^1.9)./(1+(tn./1.173474).^21.9);
Ea=En.*3.0;
E= ...
y没有定义的问题请参与 for 新手的精华帖。不够权限的话请先熟悉论坛 Dm=1;
Da=0;
end% 循环是在这结束不?
%系数矩阵A和B
循环没有结束
循环一直进行到算出y
回复 #4 youurs 的帖子
那你怎么没有end A=[ A11 A12 A13 A14 ,A21 A22 A23 A24 ,
A31 A32 A33 A34 ,
A41 A42 A43 A44 ];
粗略的看了一下，上面的那个矩阵你要表示什么意思？ 是四行，四列，还是只有16列呢？ 这样写就是16列。
第二： y没有定义可能是么有吧自己定义的函数写道equations.m中 原帖由 sigma665 于 2007-12-28 10:46 发表 http://www.chinavib.com/forum/images/common/back.gif
那你怎么没有end
哦，我没注意到，刚刚加了。
呵呵。 原帖由 re-us 于 2007-12-28 10:59 发表 http://www.chinavib.com/forum/images/common/back.gif
A=[ A11 A12 A13 A14 ,
A21 A22 A23 A24 ,
A31 A32 A33 A34 ,
A41 A42 A43 A44 ];
粗略的看了一下，上面的那个矩阵你要表示什么意思？ 是四行，四列，还是只有16列呢？ 这样写就是16列。
第二： y没有定义 ...
A应该是4*4
我改成这样了
A=[ A A A A ,
A A A A ,
A A A A ,
A A A A ];
A=-Da.*E;
A=0;
A=0;
A=Dm/0.8738;
A=Da/0.001025.*E;
A=-(Da*0.01+0.0398)/0.001025;
A=-Da/0.001025;
A=0;
A=0;
A=1/2.896;
A=-1/(0.8738.*2.896);
A=1/(0.8738.*2.896);
A=Dm*E/4;
A=0;
A=-1/(0.8738*4);
A=0;
不过还是不对
回复 #8 youurs 的帖子
什么错误，贴出来看看回复 #9 re-us 的帖子
我重新做了一下1建了一个矩阵函数
function A=mass(t,y,Da,Dm,E) %% t,y 都没有用到
A=zeros(4,4);
A(1,1)=-Da.*E;
A(1,4)=Dm/0.8738;
A(2,1)=Da/0.001025.*E;
A(2,2)=-(Da*0.01+0.0398)/0.001025;
A(2,3)=-Da/0.001025;
A(3,2)=1/2.896;
A(3,3)=-1/(0.8738.*2.896);
A(3,4)=1/(0.8738.*2.896);
A(4,1)=Dm*E/4;
A(4,3)=-1/(0.8738*4);
保存为mass.a%%为什么是*.a格式
2 微分方程
function dydt=massode(t,y,Da,Dm,E) %% 所有变量都没用到,参数没法传过来
dydt=[ A(1,1)*y(1)+A(1,4)*y(4)
A(2,1)*y(1)+A(2,2)*y(2)+A(2,3)*y(3)
A(3,2)*y(2)+A(3,3)*y(3)+A(3,4)*y(4)
A(4,1)*y(1)+A(4,3)*y(3) ];
3 主程序
t=0:0.008:2;
tn=t./(0.2+0.1555*0.8);
En=1.553174.*(tn./0.7).^1.9./(1+(tn./0.7).^1.9)./(1+(tn./1.173474).^21.9);
Ea=En.*3.0;
E=Ea+0.06;
Vl=69.36;
Qa=0;
Pa=90;
Pr=5;
if Pr>Pa
Dm=0;
Da=1;
else
Dm=1;
Da=0;
tspan=;
y0=;
options=odeset('Mass','@mass');
=ode45(@odemass,tspan,y0,options,Da,Dm,E);%% 上面是massode,这里是odemass
end
plot(t,y(:,1))
运行后
??? Undefined function or variable 'A'.
Error in ==> D:\program\matlab65\work\odemass.m
On line 2==> dydt=[ A(1,1)*y(1)+A(1,4)*y(4)
Error in ==> D:\program\matlab65\toolbox\matlab\funfun\private\odearguments.m
On line 104==> f0 = feval(ode,t0,y0,args{:});
Error in ==> D:\program\matlab65\toolbox\matlab\funfun\ode45.m
On line 155==> [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, args, ...
[ 本帖最后由 sigma665 于 2007-12-28 13:10 编辑 ]
回复 #10 youurs 的帖子
按照3楼的方法做:handshake回复 #10 youurs 的帖子
2 微分方程function dydt=massode(t,y,Da,Dm,E)
dydt=[ A(1,1)*y(1)+A(1,4)*y(4)
A(2,1)*y(1)+A(2,2)*y(2)+A(2,3)*y(3)
A(3,2)*y(2)+A(3,3)*y(3)+A(3,4)*y(4)
A(4,1)*y(1)+A(4,3)*y(3) ];
是不是也另保存了? t的变化与ode求解的时间之间有关系吗？
这个程序很费解 原帖由 youurs 于 2007-12-28 12:29 发表 http://www.chinavib.com/forum/images/common/back.gif
我重新做了一下
1建了一个矩阵函数
function A=mass(t,y,Da,Dm,E) %% t,y 都没有用到
A=zeros(4,4);
A(1,1)=-Da.*E;
A(1,4)=Dm/0.8738;
A(2,1)=Da/0.001025.*E;
A(2,2)=-(Da*0.01+0.0398)/0.001025;
...
呵呵，今天有时间，所以在罗嗦一下。需要说明的是我不知道你具体的意思，所以只能推测。
感觉你是要在3个参数变化的情况下来求解方程
你的那个if else 判断语句我也不懂，所以我只说明一种情况。
function dydt=massode(t,y)
Dm=0;
Da=1;
A=zeros(4,4);
%A(1,1)=-Da.*E;
A(1,4)=Dm/0.8738;
%A(2,1)=Da/0.001025.*E;
A(2,2)=-(Da*0.01+0.0398)/0.001025;
A(2,3)=-Da/0.001025;
A(3,2)=1/2.896;
A(3,3)=-1/(0.8738.*2.896);
A(3,4)=1/(0.8738.*2.896);
%A(4,1)=Dm*E/4;
A(4,3)=-1/(0.8738*4);
dydt=[ A(1,1)*y(1)+A(1,4)*y(4)
A(2,1)*y(1)+A(2,2)*y(2)+A(2,3)*y(3)
A(3,2)*y(2)+A(3,3)*y(3)+A(3,4)*y(4)
A(4,1)*y(1)+A(4,3)*y(3)];
存成massode.m
主程序：
for i=0:0.008:0（为了节省时间，我只循环一次）
tn=i./(0.2+0.1555*0.8);
En=1.553174.*(tn./0.7).^1.9./(1+(tn./0.7).^1.9)./(1+(tn./1.173474).^21.9);
Ea=En.*3.0;
E=Ea+0.06;
Dm=0;
Da=1;
A(1,1)=-Da.*E;
A(2,1)=Da/0.001025.*E;
A(4,1)=Dm*E/4;
tspan=;
y0=;
jsq=jsq+1;
=ode45(@massode,tspan,y0);%% 上面是massode,这里是odemass
end
plot(t,y(:,1))
最后说明，我也是新手，如有错误，原谅 原帖由 sigma665 于 2007-12-28 12:53 发表 http://www.chinavib.com/forum/images/common/back.gif
2 微分方程
function dydt=massode(t,y,Da,Dm,E)
dydt=[ A(1,1)*y(1)+A(1,4)*y(4)
A(2,1)*y(1)+A(2,2)*y(2)+A(2,3)*y(3)
A(3,2)*y(2)+A(3,3)*y(3)+A(3,4)*y(4)
A(4,1)*y(1)+A(4,3)*y(3) ];
是 ...
嗯，另存了odemass.m
[ 本帖最后由 eight 于 2007-12-29 11:24 编辑 ]
页:
[1]