声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1002|回复: 0

[编程技巧] ode计算中可以使用元胞作为变量吗?

[复制链接]
发表于 2006-7-16 14:03 | 显示全部楼层 |阅读模式

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

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

x
因为课题需要,我把变量考虑成元胞,程序如下
可调用不行,求高手指点
function myfuncircleflow1
c10=zeros(2,2);
c10(1,1)=1;
c10=c10(:);
c20=zeros(2,2);
c20=c20(:);
c30=zeros(2,2);
c30=c30(:);
c40=zeros(2,2);
c40=c40(:);
c50=zeros(2,2);
c50=c50(:);
c0={c10;c20;c30;c40;c50};
[t,y]=ode45(@myfuncircleflow,[0 8],c0);
plot(t,y);
function F=myfuncircleflow(t,c)
c1=size(4);
c2=size(4);
c3=size(4);
c4=size(4);
c5=size(4);

c={c1;c2;c3;c4;c5};


x=reshape(c1,2,2);
y=reshape(c2,2,2);
z=reshape(c3,2,2);
m=reshape(c4,2,2);
p=reshape(c5,2,2);



T=463;% K
k10=4.07e+6;%指前因子
k(1)=k10*exp(-7.88e+3/T);%速率常数
%i=1:2;
%j=1:2;

r1(1,1)=k(1)*x(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.0146);%反应速率
r1(2,1)=k(1)*x(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.0146);
r1(1,2)=k(1)*x(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.0146);
r1(2,2)=k(1)*x(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.0146);

k20=1.08e+6;
k(2)=k20*exp(-6.6e+3/T);
r2(1,1)=k(2)*y(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.000001)^0.5254;
r2(1,2)=k(2)*y(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.000001)^0.5254;
r2(2,1)=k(2)*y(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.000001)^0.5254;
r2(2,2)=k(2)*y(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.000001)^0.5254;


k30=9.8e+8;
k(3)=k30*exp(-11.16e+3/T);
r3(1,1)=k(3)*z(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.000001);
r3(1,2)=k(3)*z(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.000001);
r3(2,1)=k(3)*z(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.000001);
r3(2,2)=k(3)*z(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.000001);

k40=11.17e+9;
k(4)=k40*exp(-10.21e+3/T);
r4(1,1)=k(4)*m(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.000001);
r4(1,2)=k(4)*m(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.000001);
r4(2,1)=k(4)*m(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.000001);
r4(2,2)=k(4)*m(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.000001);

%----------------------------------------
%平衡方程
Q=0.0066;
q=0.0021;
%对PX
F=[Q*x(1,2)-Q*x(1,1)-r1(1,1);
    Q*x(1,1)-Q*x(2,1)-r1(2,1);
    Q*x(2,1)-Q*x(2,2)-r1(2,2);
    Q*x(2,2)-Q*x(1,2)-r1(1,2);
Q*y(1,2)-Q*y(1,1)+r1(1,1)-r2(1,1);
    Q*y(1,1)-Q*y(2,1)+r1(2,1)-r2(2,1);
    Q*y(2,1)-Q*y(2,2)+r1(2,2)-r2(2,2);
    Q*y(2,2)-Q*y(1,2)+r1(1,2)-r2(1,2);
Q*z(1,2)-Q*z(1,1)+r2(1,1)-r3(1,1);
    Q*z(1,1)-Q*z(2,1)+r2(2,1)-r3(2,1);
    Q*z(2,1)-Q*z(2,2)+r2(2,2)-r3(2,2);
    Q*z(2,2)-Q*z(1,2)+r2(1,2)-r3(1,2);
Q*m(1,2)-Q*m(1,1)+r3(1,1)-r4(1,1);
    Q*m(1,1)-Q*m(2,1)+r3(2,1)-r4(2,1);
    Q*m(2,1)-Q*m(2,2)+r3(2,2)-r4(2,2);
    Q*m(2,2)-Q*m(1,2)+r3(1,2)-r4(1,2);
Q*p(1,2)-Q*p(1,1)+r4(1,1);
    Q*p(1,1)-Q*p(2,1)+r4(2,1);
    Q*p(2,1)-Q*p(2,2)+r4(2,2);
    Q*p(2,2)-Q*p(1,2)+r4(1,2)];
回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-12 18:56 , Processed in 0.059154 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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