微方(ode)用解当判断条件写法

inoran

各位前辈大家好：

我有三条ＯＤＥ
我用第一条跟第二条ＯＤＥ解联立(t=0,t=0.0001)
t=0.0002开始之后就进入条件判断
做 第一条与第二条联立 ,　还有第一条与第三条联立的切换
以下是我写的 主 副程式,因为参数很多,所以可以跳过参数定义
主要是在 我多令一组未知参数矩阵 B C P Ｇ
当作我要切换的已知参数 B1 C1 P1 G1
                                       B2 C2 P2 G2

副程式是我已经加入的切换条件
程式如下:
t=ts:td:te       ts=0, td=0.0001   te=9
===========================
%解ODE副程式function dy=test2(t,y)global M K f11 e31 B33 sp h31 Q7m  z1 z2 B1 C1 P1 G1 B2 C2 P2 G2  R1 L1 Q5k B
       C P G cp  td te A

if t<0.1    V=sin(10*pi*t)if t<=0.0001

          B=B1;
          C=C1;
          P=P1;
          G=G1;

elseif  t>0.0001 && (-h31/sp).*(Q5k.*z1)*[y(1,1);y(2,1);y(3,1);y(4,1)                     %判斷條件
                                             ;y(5,1);y(6,1)]...
                                        -R1*y(17,1)-L1*y(21,1)-cp*y(13,1) > 0

             t= t-td : td : t
             B=B2;
             C=C2;
             P=P2;
             G=G2;  

elseif  t>0.0001 && (-h31/sp).*(Q5k.*z1)*[y(1,1);y(2,1);y(3,1);y(4,1)                     %判斷條件
                                             ;y(5,1);y(6,1)]...
                                        -R1*y(17,1)-L1*y(21,1)-cp*y(13,1) <= 0        

                     B=B1;
                     C=C1;
                     P=P1;
                    G=G1;
       end

else    V=0   

if  t>0.0001 && (-h31/sp).*(Q5k.*z1)*[y(1,1);y(2,1);y(3,1);y(4,1)
                                        ;y(5,1);y(6,1)]...
                                      -R1*y(17,1)-L1*y(21,1)-cp*y(13,1) > 0        

            t=t-td:td:t
             B=B2;
             C=C2;
             P=P2;
             G=G2;   

  elseif   t>0.0001 && (-h31/sp).*(Q5k.*z1)*[y(1,1);y(2,1);y(3,1);y(4,1)
                                        ;y(5,1);y(6,1)]...
                                      -R1*y(17,1)-L1*y(21,1)-cp*y(13,1) <= 0               

        B=B1;
        C=C1;
        P=P1;
        G=G1;
    end
    end

dy=[y(7:12,1)
    ;inv(M)*f11*V-inv(M)*K*y(1:6,1)...
     -inv(M)*(A.*Q7m).*(z1.*y(13,1)+z2.*y(15,1))
    ;y(17:20,1)
    ;y(21:24,1)
    ;-2.*B*y(21:24,1)-C*y(17:20,1)-(h31/sp).*P*G*y(7:12,1)]


   主程式就只是定义一堆参数 就先不看没关系,主要是这个副程式该怎么写？我看ＭＡＴＬＡＢ算ＯＤＥ是每个时间点跑一组解而我就是要在他算出来的每一组解,拿来当条件现在是我先跑前两个时间点 t=0  t=o.ooo1使用第一条ＯＤＥ与第二条ＯＤＥ联立然后t=0.0002时我就拿t=0.0001的解写个判断条件ok=>继续跑t=0.0003不ok＝＞切换成第一条与第三条联立 从t=0.0001跑到t=0.0002在跑t=0.0003  t=0.0004用t=0.0003的解再判断一次  重复这样的动作我写完是可以跑,不过感觉答案不对不知道切换会不会破坏答案的正确性？谢谢各位

[ 本帖最后由 inoran 于 2010-5-21 00:39 编辑 ]

ChaChing

看下"提问的智慧！！！！(发帖前请认真阅读)"
http://forum.vibunion.com/forum/viewthread.php?tid=21991
同待高人路过

inoran

阿！有冒犯的地方先跟您说声对不起我是想讨论看看是否有高手对於ＯＤＥ切换写法有不同的看法上面的程式是我自己写的切换程式,小弟在描述更清楚一点,就是说我有一条微方eq1想与另外第二条去解联立第二条微方因条件关系分成eq2  eq3我想检视每个时间算出来的答案去当eq2 eq3与eq1耦合的变换条件也就是说在时间开始到结束的过程在条件成立的状况下eq2 eq3会切换所以我写了上面的程式,但是使用ＯＤＥ解题器我怕在变换的过程中会有不连续产生好比说eq1 eq2在总时间长度解联立 这是连续的eq1 eq3也是一样但是在整个时域的每个时间点如果有  eq2 与eq3切换跟eq1耦合去联立解的时候我想法是会有不连续的现象,所以用ＯＤＥ解题器,是否就无法实现我这样的作法或者说其实可以做到？,我就ＰＯ了我写的程式想跟大家讨论看看如果有不礼貌或是发言不慎的地方,请见谅,真的很对不起＠＠

ChaChing

别误会LZ那有冒犯!?
可能个人水平专业有限, 感觉LZ的问题有点乱, 说真的还没空认真看! :@L
所以就仅是善意提醒, 别介意!:loveliness:

messenger

建议你还是自己编程解决吧，自己编程也不是很难。
Matlab因为采用自适应步长，实现你这个要求很难。就好像一个全自动照像机，要实现手动调整一样。