关于四阶五级Runge-Kutta变步长解微分方程的问题

youurs

这是我的程序
%算变量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,Qa,Pa,Pr]')
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=[B1 B2 B3 B4];
dydt=A*y+B;
tspan=[0 2];
y0=[69.36;0;90;5];
options=odeset('Vl','Qa','Pa','Pr');
[t,y]=ode45(@equations,tspan,y0,options);
plot(t,y(:,1))

为什么运行后说y是没有定义的？
ps：options=odeset('Vl','Qa','Pa','Pr')这个设置对吗？？  谢谢！！

[ 本帖最后由 ChaChing 于 2009-4-17 08:41 编辑 ]

sigma665

??? Attempt to execute SCRIPT function as a function:
D:\Program Files\MATLAB\R2007b\toolbox\matlab\lang\function.m

在红色加了end后
运行错误如上

Function dydt=equations(t,y)有错误

eight

原帖由 youurs 于 2007-12-28 08:48 发表

                               
登录/注册后可看大图

这是我的程序
%算变量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 新手的精华帖。不够权限的话请先熟悉论坛

youurs

Dm=1;
Da=0;
end% 循环是在这结束不?
%系数矩阵A和B

循环没有结束
循环一直进行到算出y

sigma665

那你怎么没有end

re-us

A=[ A11 A12 A13 A14 ,
A21 A22 A23 A24 ,
A31 A32 A33 A34 ,
A41 A42 A43 A44 ];
粗略的看了一下，上面的那个矩阵你要表示什么意思？ 是四行，四列，还是只有16列呢？ 这样写就是16列。
第二： y没有定义可能是么有吧自己定义的函数写道equations.m中

youurs

原帖由 sigma665 于 2007-12-28 10:46 发表

                               
登录/注册后可看大图

那你怎么没有end

哦，我没注意到，刚刚加了。
呵呵。

youurs

原帖由 re-us 于 2007-12-28 10:59 发表

                               
登录/注册后可看大图

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[1,1] A[1,2] A[1,3] A[1,4] ,
A[2,1] A[2,2] A[2,3] A[2,4] ,
A[3,1] A[3,2] A[3,3] A[3,4] ,
A[4,1] A[4,2] A[4,3] A[4,4] ];

A[1,1]=-Da.*E;
A[1,2]=0;
A[1,3]=0;
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[2,4]=0;
A[3,1]=0;
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,2]=0;
A[4,3]=-1/(0.8738*4);
A[4,4]=0;

不过还是不对

re-us

什么错误，贴出来看看

youurs

我重新做了一下
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=[0 2];
y0=[69.36;0;90;5];
options=odeset('Mass','@mass');
[t,y]=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 编辑 ]

花如月

按照3楼的方法做:handshake

sigma665

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) ];
是不是也另保存了?

无水1324

t的变化与ode求解的时间之间有关系吗？
这个程序很费解

re-us

原帖由 youurs 于 2007-12-28 12:29 发表

                               
登录/注册后可看大图

我重新做了一下
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=[0 2];
y0=[69.36;0;90;5];
jsq=jsq+1;
[t,y]=ode45(@massode,tspan,y0);%% 上面是massode,这里是odemass
end
plot(t,y(:,1))

最后说明，我也是新手，如有错误，原谅

youurs

原帖由 sigma665 于 2007-12-28 12:53 发表

                               
登录/注册后可看大图

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 编辑 ]