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

cichishia

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛

! u8 w4 ^/ _, f; h. B# L- S这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
case 'creep & recovery'
        stress0=arg(1);
: z3 W' W) R/ Y, V! b6 r( v        risetime=arg(2);
        creeptime=arg(3);
        if t<risetime
8 c4 _7 `- w) s3 \  u            stressrate=stress0/risetime;
3 j+ H5 b* p& \) b            stress=stressrate*t;
        elseif t<creeptime
1 }' Z, h; a& A& ]) q$ n            stressrate=0;
            stress=stress0;
        elseif t<creeptime+risetime
            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
. |; w# i- _; R  Z4 c- ^; m        else
            stressrate=0;
3 |5 H3 x- V) ~9 X- o( F6 S            stress=0;
# k9 f7 \/ m* {& w0 `4 z        end
+ _# g* b+ P/ a  o* H  m        gam1=y(1);
$ g- ~. f. k/ r. B. K5 |+ ?( R$ B        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);

  }, g& o! R/ A. E7 ?1 q

" K) s  W0 ^  L4 u" g' t* U. Rclear all
clc
warning off
) H: Y5 S0 p/ j+ J
6 l+ j& m% F0 u8 y1 ueta1=10;  % [Pa.s]
  z) q8 T" U5 Q' E, U2 v6 veta2=20;  % Pa.s]
! c( G: _% z) q' @* E) QG1=2;   % [Pa]
+ R2 O! X* X( D6 {$ oG2=5;   %[Pa]
( F8 q9 A0 q4 s" v; ]yieldstress=5;  % [Pa]
) J" x/ h: M) a' ?: ?contactstrain=5;

+ H2 a5 _" y- K5 v, w# j" D  w8 I
3 I7 A$ r2 O9 Q) G+ D  sflow='creep & recovery';
* t4 W9 T" W/ m9 J5 O& ustress0=1:2:10;
risetime=0.01;
creeptime=100;
  ?: }( r8 V, Z: K) s6 mfor i=1:length(stress0)
    [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
5 J$ ?! l2 O3 j    a{i}=strcatnum2str(stress0(i),'(Pa)');
# a, M2 o& m$ R* _0 @; a$ Aend
/ [& l: c* T9 G/ V( p/ Qxlin([1e-4 1000]);   
* E! K- J. K8 I& J' _6 I7 z4 G* Pxlabel('Time(s)');      ylabel('Shear strain(-)');
% {! }: @4 A9 z8 q" V; P+ rlegend(a)


: `- ~- n( C6 Q6 I7 [, M- l& `4 F
6 I% V* r4 k- k: q# |  o错误使用 odearguments (第 93 行)
! i9 ?; n9 B: n( l7 n5 Q1 [MYWORK 必须返回列向量。

1 v* @  l% D* z9 o$ i. e出错 ode15s (第 152 行)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

出错 Untitled2 (第 20 行)
, D0 |+ R# X% e8 e* u    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
7 K& T& J  I6 s
. e6 G* M- p9 l. [
$ n4 D/ d8 Y8 i4 C& _( j+ D# _感谢大神们！
- ^( I+ X0 r+ l1 ~* m, v, d

llbnmo

clear all
, I9 U; O- v: lclc
warning off
4 @3 c! B0 |# x
7 l9 K1 q5 t0 g! reta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
$ V$ @! V, w' [* c4 m! ]G1=2;   % [Pa]
G2=5;   %[Pa]
yieldstress=5;  % [Pa]
contactstrain=5;
- Y7 B! S4 u9 Y+ F& f. s
- r! _0 U  y9 `
. P) {# A. j& f/ V: U8 A+ Z" Eflow='creep & recovery';
5 Q" L6 P6 H4 p0 X4 {, n4 Ystress0=1:2:10;
risetime=0.01;
creeptime=100;
1 p- f; L2 z& A( S. Afor i=1:length(stress0)
! _8 ^2 V7 W! ^1 m& d5 X    [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
4 i! m* e2 u2 P2 w' X: y4 ?9 A$ B    a{i}=strcat(num2str(stress0(i)),'(Pa)');
+ U: `5 i8 \( F- I7 C5 H$ o% x; Nend
xlim([1e-4 1000]);
2 @0 J6 q0 ?8 ^. ?. V+ H0 qxlabel('Time(s)');      ylabel('Shear strain(-)');
1 o' k& [% x/ @& l7 Blegend(a)
! w( f# o0 i+ g0 g+ @function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
    case 'creep & recovery'
        stress0=arg(1);
        risetime=arg(2);
        creeptime=arg(3);
+ y; H. ]* c5 M        if t<risetime
            stressrate=stress0/risetime;
) n% p: O  r, {# q6 k9 h" m            stress=stressrate*t;
  f9 H# j& t$ h" m9 w! i0 J        elseif t<creeptime
& P8 G$ B& G( \5 |- @9 Y            stressrate=0;
& x/ J5 y+ V& w: ^$ _; j            stress=stress0;
+ |, ?' j7 u* w$ C" j        elseif t<creeptime+risetime
            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
            stressrate=0;
            stress=0;
        end
        gam1=y(1);
        gam2=y(2);
        s(1)=gam2;
% M4 j+ w# T7 s" I( A        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);
' S7 ~+ q0 }1 s! T( L) ~        s=s';
end
8 z2 j& t' ~' e& h- _  q, Qend

小白的白

楼上正解！

IBB-EUT

来学习一下