matlab解常微分一直报错返回列向量

cichishia

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
$ I6 F. Y! }' e# O: p" W
这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
5 K8 X+ o5 n* L( s( Vswitch flow
case 'creep & recovery'
+ [( s6 {* `- C: ]& }4 x        stress0=arg(1);
        risetime=arg(2);
* d/ {" G/ a- O; A        creeptime=arg(3);
* D- ~" R9 k5 p  W        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
            stress=stress0;
7 V7 S; I3 v, H+ ?% L        elseif t<creeptime+risetime
4 ^( B% o) [/ Q8 Y5 c$ \$ h& b            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
2 V% K( e2 m' {- c1 a. E* {, X        else
- t. s; @' }5 q4 }1 `, R( l3 \            stressrate=0;
            stress=0;
' d! S  z" {4 r        end
; c0 Q" j7 h' S0 F) X8 D        gam1=y(1);
        gam2=y(2);
        s(1)=gam2;
        s(2)=G1/(eta1*eta2*((abs(gam1)-contactstrain)>0))*((stress-eta2*gam2*((abs(gam1)-contactstrain)>0)-sign(stress-eta2*gam2*((abs(gam1)-contactstrain)>0))*yieldstress)*((abs(stress-eta2*gam2*((abs(gam1)-contactstrain)>0))-yieldstress)>0)-eta1*gam2+eta1/G1*stressrate);
. F0 n! v1 s2 N1 e! W- v# O3 _/ P1 A, E
9 x# v5 e. u! A4 J) e' g, p
7 [3 K4 Q0 j" y% v3 `3 C
clear all
& v3 t6 n2 P! C1 Mclc
1 a2 Q$ O; ?; Y. c6 G; o% N0 u# Awarning off

5 |% P* |- f/ C: veta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
& f% }  ^7 N/ D8 mG1=2;   % [Pa]
  \) y5 U5 t4 G: C9 j; l$ tG2=5;   %[Pa]
/ H9 P9 C) E- |. t6 i( c& F" ^5 \yieldstress=5;  % [Pa]
contactstrain=5;


5 F9 f1 t) v9 qflow='creep & recovery';
0 j+ @0 l, W- vstress0=1:2:10;
; D' F& y5 h2 p7 m5 {% V3 ~) z' v6 Erisetime=0.01;
creeptime=100;
+ Y9 q! ^+ z3 _+ Z* l9 d/ ]for i=1:length(stress0)
( z2 [9 _* B/ q0 A. C4 L6 V5 E- p% o    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
! P3 P5 o6 v" c    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcatnum2str(stress0(i),'(Pa)');
( C: D% |7 D- a* n' s1 }5 Kend
xlin([1e-4 1000]);   
xlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)

$ v& O# T' T, w, Y/ |1 ^

& r7 s0 Z& ?; H错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。
7 H! _1 T# J: l! t& a6 h9 ?
出错 ode15s (第 152 行)
7 q3 {7 w/ X6 X    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

出错 Untitled2 (第 20 行)
" I# s6 I6 l3 T- o! n* j# D3 T/ ]    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);


感谢大神们！
+ `$ A, g# _- ^! D  r# b

llbnmo

clear all
1 T0 ]4 H, X5 K' d# T% m" lclc
, P' w& {8 M" u5 V& _warning off

& ~1 L% H& C$ U0 O/ n4 A0 [, zeta1=10;  % [Pa.s]
) [9 |( i/ d+ Aeta2=20;  % Pa.s]
G1=2;   % [Pa]
) V4 m! H  o! o! M( S6 G2 `, lG2=5;   %[Pa]
yieldstress=5;  % [Pa]
: `1 C" s9 T1 Y. @: L& tcontactstrain=5;


& s  N) O2 _7 u& e, Mflow='creep & recovery';
/ ^( t  V( y' K" Xstress0=1:2:10;
' r0 Z# {7 |4 X5 Irisetime=0.01;
% E; T+ ?7 c, a$ @" W. n# N# qcreeptime=100;
5 B5 x( w- V; x5 dfor i=1:length(stress0)
5 i. o( l4 Q% D& ?. I. c8 W    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
: m3 ^1 j+ Q! R! ]2 B    a{i}=strcat(num2str(stress0(i)),'(Pa)');
end
" l0 Q5 |1 [$ E, u2 C/ X, Nxlim([1e-4 1000]);
xlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
  v& ]9 u6 `' t7 q( m$ b$ cswitch flow
$ ~& w5 Y) p1 V$ H# Y0 H    case 'creep & recovery'
9 K! S! D% j. Y        stress0=arg(1);
        risetime=arg(2);
4 [# s1 {6 X' O  ]6 z; I9 m        creeptime=arg(3);
7 v4 T! V6 ~6 J/ @2 e0 p        if t<risetime
9 X; O  {4 r) B( }            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
) ]6 m9 |- O2 B* M; D8 v            stressrate=0;
% r( A7 D2 B3 l. W/ n# @            stress=stress0;
! y- X  n- ?# k        elseif t<creeptime+risetime
            stressrate=stress0/risetime;
( n3 q8 B+ ?& w; F8 t            stress=stress0-stressrate*(t-creeptime);
( k- p3 ?0 A  l" F/ v        else
- D0 z3 ^) l: R5 |: |( S            stressrate=0;
/ R, S4 V0 c  X, m5 I) `0 Z            stress=0;
        end
        gam1=y(1);
1 U$ c, x+ E2 m0 @4 \1 C( T        gam2=y(2);
        s(1)=gam2;
- t/ H$ U) y! @* Z  d+ w% F& x        s(2)=G1/(eta1*eta2*((abs(gam1)-contactstrain)>0))*((stress-eta2*gam2*((abs(gam1)-contactstrain)>0)-sign(stress-eta2*gam2*((abs(gam1)-contactstrain)>0))*yieldstress)*((abs(stress-eta2*gam2*((abs(gam1)-contactstrain)>0))-yieldstress)>0)-eta1*gam2+eta1/G1*stressrate);
        s=s';
end
: R$ t0 D4 h$ ^6 E! T, [/ Mend

小白的白

楼上正解！

IBB-EUT

来学习一下