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

cichishia

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

这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
: [* r* ?' Q! l# o6 s: j7 p0 U0 W* f9 \switch flow
# y5 _1 Y$ Q  L) J5 zcase 'creep & recovery'
        stress0=arg(1);
        risetime=arg(2);
4 L  N+ j8 l& i/ [  e        creeptime=arg(3);
8 ?( m, [, F; O. U        if t<risetime
/ C7 E. k" y* n8 R            stressrate=stress0/risetime;
1 t) _$ C% P/ f- q5 ^. ^            stress=stressrate*t;
' X  d% Q7 k/ g2 M" c- N4 u" V9 m        elseif t<creeptime
            stressrate=0;
            stress=stress0;
7 E  H! V% ?; ]/ y- Z* C+ z        elseif t<creeptime+risetime
- U7 F# j0 i  B            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
) d+ l; `( D$ z        else
            stressrate=0;
            stress=0;
        end
        gam1=y(1);
        gam2=y(2);
        s(1)=gam2;
$ F! [% Z8 l9 S/ z        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);
. C2 }( p6 u8 Q" r  m
: G. d, Y% b( ~9 P, L6 |+ n4 m! s

( O; X, C3 b) k/ O( b# oclear all
. Z0 f' x9 z6 B+ K  `+ z/ V/ Nclc
8 j' _  U# H0 Awarning off
  X9 [* \+ b; |# p
eta1=10;  % [Pa.s]
5 N, v0 x6 `3 O( R7 x3 eeta2=20;  % Pa.s]
G1=2;   % [Pa]
G2=5;   %[Pa]
, i) C& Q  C8 K6 qyieldstress=5;  % [Pa]
. {- P$ ?5 p1 j, ?2 D( ^4 |7 {" jcontactstrain=5;
4 b* W6 I" k$ |" J3 s7 k; a6 K

flow='creep & recovery';
stress0=1:2:10;
risetime=0.01;
& g$ J% {" ?% U3 R. _6 v! m" U' rcreeptime=100;
for i=1:length(stress0)
# }1 v+ [/ {" W" Y" h. n7 Q    [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
, H. o8 |3 h! X) x    a{i}=strcatnum2str(stress0(i),'(Pa)');
6 s3 D7 }3 `# H* H9 I) iend
1 Y# F* _( C# M9 I2 f" B4 p" M& txlin([1e-4 1000]);   
xlabel('Time(s)');      ylabel('Shear strain(-)');
( F9 G& _! [9 v0 s6 _" v# ]( l3 Zlegend(a)
4 _3 y) n; z( E5 G- v/ y8 [
% Y- R, p& Z. D7 ?$ [) P: z$ `

错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。
; l+ d, S9 [) {4 u/ w; ~
出错 ode15s (第 152 行)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
' p0 {: O* w0 o* Y8 h; r
% A0 @) H! H* y; D出错 Untitled2 (第 20 行)
2 @! n; ~* G# H- Y$ p    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);


  J0 c  l. a2 g& ^+ g) M感谢大神们！
2 v$ Q5 ]- h$ b4 D8 C7 r: y

llbnmo

clear all
: M& }! E2 y' M) v2 k- f+ Rclc
warning off
5 k2 I, z7 t. L9 d: Q0 P9 z
eta1=10;  % [Pa.s]
/ w5 L: }: J1 E+ U( C% ?& H" ^6 Leta2=20;  % Pa.s]
/ p5 d* B- Y' A7 q/ `- H- n; W+ L8 iG1=2;   % [Pa]
/ q7 e/ z4 W2 C8 TG2=5;   %[Pa]
yieldstress=5;  % [Pa]
. U9 ?7 x! c9 I' hcontactstrain=5;

; v: V1 l. `- a3 B. @& _
flow='creep & recovery';
5 }/ c3 H! J, C+ b9 @4 Xstress0=1:2:10;
risetime=0.01;
, c) W0 p5 U) a* {% Screeptime=100;
. x5 d7 U" V4 z. t% {for i=1:length(stress0)
; v; V3 G; A) w# V9 M# n& m/ A    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
8 J, Y9 f4 f# L5 k5 R    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcat(num2str(stress0(i)),'(Pa)');
end
- X; }+ `: f+ d( t5 r1 [! dxlim([1e-4 1000]);
& c; B  B' W. t" yxlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
: q  a5 d) ^( {( E4 Efunction s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
9 T+ Z5 m0 ]4 V% Y! Iswitch flow
# R3 t' \5 W" L3 z8 y    case 'creep & recovery'
0 O) s. I: r+ j8 W: r$ G8 A" \        stress0=arg(1);
# p/ S+ F0 ?: w( k8 Y1 f        risetime=arg(2);
& ?8 v, O, x+ v/ R/ _        creeptime=arg(3);
        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
  q0 d" \3 x7 e! x- o; F) z; ]        elseif t<creeptime
            stressrate=0;
7 t1 ~  g7 b# b/ n7 B            stress=stress0;
% Y- P5 O- y2 W9 v7 D2 ]        elseif t<creeptime+risetime
            stressrate=stress0/risetime;
+ X# {. s: h2 n! h; Z& E- k. N            stress=stress0-stressrate*(t-creeptime);
        else
' d' N) i5 ?8 H; ^; B9 G' W% U            stressrate=0;
4 _1 x* B/ c, t9 X+ H, q3 O            stress=0;
        end
        gam1=y(1);
2 D( d7 K, u  r" b2 A4 S/ c        gam2=y(2);
/ h, }) n" Z/ X6 P! D% {; P        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);
        s=s';
end
end

小白的白

楼上正解！

IBB-EUT

来学习一下