应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » matlab不报错但结果不对啊，
61/1返回列表
查看: 1067  |  回复: 5
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

liutinglv

金虫 (著名写手)

[求助] matlab不报错但结果不对啊，

同胞们，小硕一事不明，就是matalab在求解微分方程的时候不报错但是，运行结果不对啊和想像中的很大偏差
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼2012-03-28 16:24:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ldoop

铁杆木虫 (著名写手)

【答案】应助回帖

感谢参与，应助指数 +1
建议把code贴上来，检查一下里面是否表达式写错了之类的。
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

liutinglv
2楼2012-03-28 17:16:47
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

liutinglv

金虫 (著名写手)
送鲜花一朵
dbb627: 编辑内容 2012-03-28 21:56
引用回帖:
2楼: Originally posted by ldoop at 2012-03-28 17:16:47:
建议把code贴上来，检查一下里面是否表达式写错了之类的。

CODE:
function gudingchuang6
clear all
clc
global A B C E F H I J p R e h  X Y V y A1 B0  W0 B1 M N a b A2 G R0 R1 J1 Cb3

A=[7/3  -8/3  1/3
  -8/3  16/3 -8/3
   1/3  -8/3  7/3];
B=[-1/2  2/3  -1/6
   -2/3   0    2/3
    1/6 -2/3  1/2];
C=[2/15 1/15 -1/30
   1/15 8/15 1/15
  -1/30 1/15 2/15];
F=[24.00 -37.18 25.18 -12.00
   16.39 -24.00 12.00 -4.392
   -4.39 12.00 -24.00 16.39
  -12.00 25.18 -37.18 24.00];
E=[-7.0000  8.1960  -2.1960  -1.0000
    -2.7320 1.7320 1.7320   -0.7320
    0.7320 -1.7320 -1.7320  2.7320
    -1.0000 2.1960 -8.1960  7.0000];
G=[0.002983, 0.007894, 0.01818, 0.03642, 0.06350, 0.09635, 0.1272, 0.147, 0.147, 0.1272, 0.09635, 0.06350, 0.03642, 0.01818, 0.007894, 0.002989];
H=10*[0.006567, 0.005264, 0.004392, 0.003768, 0.0033, 0.002935, 0.002643, 0.002403, 0.002204, 0.002035, 0.00189, 0.001764, 0.001654, 0.001557, 0.001471, 0.001394];
I=[3.2948, 2.8103, 2.4817, 2.2308, 2.0372, 1.8789, 1.7536, 1.6462, 1.5548, 1.4756, 1.4066, 1.3449, 1.2902, 1.2409, 1.197, 1.1569];
J=[20.44, 21.75, 23.02, 24.12, 25.15, 26.08, 27.03, 27.91, 28.74, 29.54, 30.32, 31.06, 31.78, 32.47, 33.15, 33.81];
R=[0.0001, 0.1057, 0.3943, 0.5, 0.6057, 0.8943, 1] ;
a=39;
b=81;
p=16.31;
e=0.6;
h=1/2;
X=(1/(p*h))*A;
y=sum(G.*I);
Y=h*y*C;
U=X+B+Y;
V=(-h)*C;
A1=sparse(5,5);
B0=sparse(5,5);
A1(1,1:3)=U(1,1:3);
B0(1,1:3)=V(1,1:3);
for i=2:2:4,
    A1(i,(i-1):(i+1))=U(2,1:3);
    B0(i,(i-1):(i+1))=V(2,1:3);
end
for i=3;
    A1(i,(i-2):(i-1))=U(3,1:2);
    A1(i,i)=U(3,3)+U(1,1);
    A1(i,(i+1):(i+2))=U(3,2:3);
    B0(i,(i-2):(i-1))=V(3,1:2);
    B0(i,i)=V(3,3)+V(1,1);
    B0(i,(i+1):(i+2))=V(3,2:3);
end
    A1(5,3:5)=U(3,1:3);
    B0(5,3:5)=V(3,1:3);
    B1=inv(B0);
A2=B1*A1;
M=sparse(7,7);
R0=2*ones(1,7);
R1=R0./R;
M(1,1:4)=E(1,1:4);
for i=2:3:5;
M(i,(i-1):(i+2))=2*R1(i)*E(2,1:4)+4*F(2,1:4);
end;
for i=3:3:6;
M(i,(i-2):(i+1))=2*R1(i)*E(3,1:4)+4*F(3,1:4);
end;
for i=4;
M(4,1)=E(4,1);
M(4,2)=E(4,2);
M(4,3)=E(4,3);
M(4,4)=E(4,4)-E(1,1);
M(4,5)=-E(1,1);
M(4,6)=-E(1,2);
M(4,7)=-E(1,1);
end
M(7,4:7)=E(4,1:4);
W0=zeros(1,325); ;
[t,W]=ode45(@Equations,[0,2],W0);

Cp(:,1:320)=W(:,1:320);
Cb(:,1:5)=W(:,321:325);
  for i=1:80,
    for j=1:length(t),
        Cp4(j,i)=((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp(j,i))+(-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp(j,(80+i)))+ ...
  (M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp(j,(160+i))+(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4))*Cp(j,240+i))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))...
       +M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1);
    end
end
for i=1:80,
   for j=1:length(t),
       Cp1(j,i)=-(M(1,2)*Cp(j,i)+M(1,2)*Cp(j,(80+i)))/M(1,1)-(M(1,4)/M(1,1))*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp(j,i)+...
           (-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp(j,(i+80))+(M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp(j,(i+160))+...
(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4)))*Cp(j,i+240)))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))+M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1));
   end
end
for j=1:length(t),
for i=1:5:76,
    Cb3(j,i)=Cb(j,1);
end
end
for j=1:length(t),
for i=2:5:77,
    Cb3(j,i)=Cb(j,2);
end
end
for j=1:length(t),
for i=3:5:78,
    Cb3(j,i)=Cb(j,3);
end
end
for j=1:length(t),
for i=4:5:79,
    Cb3(j,i)=Cb(j,4);
end
end
for j=1:length(t),
for i=5:5:80,
    Cb3(j,i)=Cb(j,5);
end
end
for j=1:length(t),
     for m=1:16,
          for i=1:5,
           J1(j,(5*(m-1)+i))=J(m);
          end
     end
end
for i=1:80,
      for j=1:length(t),
        Cp7(j,i)=(J1(j,i)*Cb3(j,i)-M(7,5)*Cp(j,(i+160))-M(7,6)*Cp(j,(i+240))/(J1(j,i)+M(7,7))-M(7,4)*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp(j,i))...
    +(-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp(j,(i+80))+(M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp(j,(i+160))+(M(4,6)/M(4,4)-M(4,7)*M(7,6) ...
    /(M(7,7)*M(4,4)))*Cp(j,(i+240))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4)))+M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1)))/(J1(j,i)+M(7,7));
      end
end
Cp=[Cp Cp1 Cp4 Cp7];
disp([W t]);
function [dWdt]=Equations(t,W)
global H I J p R e h  B1 M N a b A2 G  f g J1  N11 Cb3
Cp(1:320)=W(1:320);
Cb(1:5)=W(321:325);
for i=1:5:76,
    Cb3(i)=Cb(1);
end
for i=2:5:77,
    Cb3(i)=Cb(2);
end
for i=3:5:78,
    Cb3(i)=Cb(3);
end
for i=4:5:79,
    Cb3(i)=Cb(4);
end
for i=5:5:80,
    Cb3(i)=Cb(5);
end
for j=1:16,
          for i=1:5,
           J1(5*(j-1)+i)=J(j);
          end
     end
for i=1:80,
Cp7(i)=(J1(i)*Cb3(i)-M(7,5)*Cp((i+160))-M(7,6)*Cp((i+240))/(J1(i)+M(7,7))-M(7,4)*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp(i))...
    +(-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp((i+80))+(M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp((i+160))+(M(4,6)/M(4,4)-M(4,7)*M(7,6) ...
    /(M(7,7)*M(4,4)))*Cp((i+240))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4)))+M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1)))/(J1(i)+M(7,7));
end
for m=1:16,
    for i=1:5,
         N11(m,i)=M(2,2)*Cp(5*(m-1)+i)+M(2,3)*Cp((5*(m-1)+80+i))+M(2,1)*(-(M(1,2)*Cp((5*(m-1)+i))+M(1,2)*Cp((80+(5*(m-1)+i)))/M(1,1)-(M(1,4)/M(1,1))*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp((5*(m-1)+i))+...
           (-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp(((5*(m-1)+i)+80))+(M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp(((5*(m-1)+i)+160))+...
(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4)))*Cp((5*(m-1)+i)+240)))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))+M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1)))+...
   M(2,4)*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp((5*(m-1)+i)))+(-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp((80+(5*(m-1)+i))))+ ...
  (M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp((160+(5*(m-1)+i)))+(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4))*Cp(240+(5*(m-1)+i)))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))...
       +M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1)));
    end
end
for i=1:5
    for m=1:16,
        N11(m,(i+5))=M(3,2)*Cp(5*(m-1)+i)+M(3,3)*Cp((5*(m-1)+80+i))+M(3,1)*(-(M(1,2)*Cp((5*(m-1)+i))+M(1,2)*Cp((80+(5*(m-1)+i)))/M(1,1)-(M(1,4)/M(1,1))*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp((5*(m-1)+i))+...
           (-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp(((5*(m-1)+i)+80))+(M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp(((5*(m-1)+i)+160))+...
(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4)))*Cp((5*(m-1)+i)+240)))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))+M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1)))+...
   M(3,4)*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp((5*(m-1)+i)))+(-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp((80+(5*(m-1)+i))))+ ...
  (M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp((160+(5*(m-1)+i)))+(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4))*Cp(240+(5*(m-1)+i)))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))...
       +M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1)));
    end
end


for i=1:5
    for m=1:16,
        N11(m,(i+10))=M(5,5)*Cp((5*(m-1)+160+i))+M(5,6)*Cp((5*(m-1)+240+i))+M(5,4)*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp((5*(m-1)+i)))+(-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp((80+(5*(m-1)+i))))+ ...
  (M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp((160+(5*(m-1)+i)))+(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4))*Cp(240+(5*(m-1)+i)))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))...
       +M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1))+M(5,7)*Cp7(5*(m-1)+i);
    end
end
for i=1:5
    for m=1:16,
        N11(m,(i+15))=M(6,5)*Cp((5*(m-1)+160+i))+M(6,6)*Cp((5*(m-1)+240+i))+M(6,4)*(((-M(4,1)*M(1,2)/(M(1,1)*M(4,4))+M(4,2)/M(4,4))*Cp((5*(m-1)+i)))+(-M(4,1)*M(1,3)/((M(1,1)*M(4,4))+M(4,3)/M(4,4))*Cp((80+(5*(m-1)+i))))+ ...
  (M(4,5)/M(4,4)-M(4,7)*M(7,5)/(M(7,7)*M(4,4)))*Cp((160+(5*(m-1)+i)))+(M(4,6)/M(4,4)-M(4,7)*M(7,6)/(M(7,7)*M(4,4))*Cp(240+(5*(m-1)+i)))/(M(4,1)*M(1,4)/(M(1,1)*M(4,4))...
       +M(4,7)*M(7,4)/(M(7,7)*M(4,4))-1))+M(6,7)*Cp7(5*(m-1)+i);
    end
end

Cp8=reshape(Cp7,5,16);
H1=G.*I;
dCbdt(1)=A2(1,1:5)*Cb'+B1(1,1)*Cb(1)-B1(1,1)+sum(H1.*Cp8(1,:));
for k=2:5,
dCbdt(k)=A2(k,1:5)*Cb'+sum(H1.*Cp8(k,:));
end
H2(1:16)=G.*H;
for i=1:16,
     for k=1:5,
        for m=1:4,
          dCpdt(80*(m-1)+5*(i-1)+k)=1/((a*(1-e)/(1+b*Cp(80*(m-1)+5*(i-1)+k))^2)+e)*(H2*N11(:,(5*(m-1)+k)));
        end
     end
end
  dWdt=[dCpdt dCbdt]';

辛苦高手帮忙看看

[ Last edited by dbb627 on 2012-3-28 at 21:56 ]
赞一下 回复此楼
3楼2012-03-28 20:59:08
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

rex881026

木虫 (小有名气)
表示膜拜
回复此楼
4楼2012-03-28 21:59:39
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

joy4210

金虫 (小有名气)

【答案】应助回帖

★
感谢参与，应助指数 +1
zhangguangping: 金币+1, 谢谢建议！ 2012-03-31 00:10:45
那看你的表达式与程序中一致么？
赞一下 回复此楼
5楼2012-03-30 08:22:13
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

695710829

新虫 (初入文坛)

【答案】应助回帖

感谢参与，应助指数 +1
你的公式比较长，你确定没编错？还有建议你检查检查循环部分哦。
赞一下 回复此楼
6楼2012-03-31 12:50:50
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 liutinglv 的主题更新
61/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭