声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3142|回复: 14

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

[复制链接]
发表于 2007-12-28 08:48 | 显示全部楼层 |阅读模式

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

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

x
这是我的程序
%算变量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 编辑 ]
回复
分享到:

使用道具 举报

发表于 2007-12-28 09:27 | 显示全部楼层

回复 #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)有错误
发表于 2007-12-28 09:39 | 显示全部楼层
原帖由 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 新手的精华帖。不够权限的话请先熟悉论坛
 楼主| 发表于 2007-12-28 10:42 | 显示全部楼层
Dm=1;
Da=0;
end% 循环是在这结束不?
%系数矩阵A和B

循环没有结束
循环一直进行到算出y
发表于 2007-12-28 10:46 | 显示全部楼层

回复 #4 youurs 的帖子

那你怎么没有end
发表于 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没有定义可能是么有吧自己定义的函数写道equations.m中
 楼主| 发表于 2007-12-28 11:10 | 显示全部楼层
原帖由 sigma665 于 2007-12-28 10:46 发表
那你怎么没有end


哦,我没注意到,刚刚加了。
呵呵。
 楼主| 发表于 2007-12-28 11:13 | 显示全部楼层
原帖由 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;

不过还是不对
发表于 2007-12-28 11:21 | 显示全部楼层

回复 #8 youurs 的帖子

什么错误,贴出来看看
 楼主| 发表于 2007-12-28 12:29 | 显示全部楼层

回复 #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=[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 编辑 ]
发表于 2007-12-28 12:49 | 显示全部楼层

回复 #10 youurs 的帖子

按照3楼的方法做:handshake
发表于 2007-12-28 12:53 | 显示全部楼层

回复 #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) ];
是不是也另保存了?
发表于 2007-12-28 13:44 | 显示全部楼层
t的变化与ode求解的时间之间有关系吗?
这个程序很费解

评分

1

查看全部评分

发表于 2007-12-28 15:34 | 显示全部楼层
原帖由 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))

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

评分

1

查看全部评分

 楼主| 发表于 2007-12-29 11:18 | 显示全部楼层
原帖由 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 编辑 ]
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-28 12:20 , Processed in 0.076991 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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