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

cichishia

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
4 H8 i+ N+ w$ W! x
& O6 T( ^& h, A( r$ Q$ `$ g: q/ K这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
) d9 K% `! T, k- s: p9 v# G0 {case 'creep & recovery'
        stress0=arg(1);
6 J. t2 K, h$ f' p        risetime=arg(2);
        creeptime=arg(3);
        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
  P6 S; v. N, X, m6 ?        elseif t<creeptime
            stressrate=0;
8 p) ^; o( l6 R7 w; ^' K            stress=stress0;
        elseif t<creeptime+risetime
# y* X  g4 w$ t0 O' a4 |2 S8 _            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
& L+ Z& a4 ?9 c            stressrate=0;
. D/ {; i7 _" l            stress=0;
$ e: v# U' @( v- i6 h- G        end
        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);
1 U& q: t, ?- `4 {) X+ K, |


clear all
' M* U& N% |8 Q, @7 mclc
warning off
( a" [( n" W( Q3 |- }
eta1=10;  % [Pa.s]
1 h0 ~$ z5 r; @4 e' ~eta2=20;  % Pa.s]
6 L. N$ H- h5 X5 bG1=2;   % [Pa]
G2=5;   %[Pa]
yieldstress=5;  % [Pa]
" h' k0 n; U: r% Ncontactstrain=5;
2 L$ @; }8 J; v5 k7 I

flow='creep & recovery';
stress0=1:2:10;
% ?( }. _5 |# C" u( ?9 orisetime=0.01;
; X; A  `% V% b+ f# [creeptime=100;
for i=1:length(stress0)
, j7 J; S3 g# {4 \8 e    [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
    a{i}=strcatnum2str(stress0(i),'(Pa)');
9 V# G/ e4 V  w4 f) Wend
4 N2 ~3 b, X! q' [/ M( c1 R- N5 uxlin([1e-4 1000]);   
0 W# o: f, i1 s6 O) ~  j  \4 `5 nxlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
: L" F2 ~6 r9 U

  x$ M6 H: [3 z! N7 @7 \
错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。
9 I& ]( x7 p3 f0 m- y6 H) [
出错 ode15s (第 152 行)
* `; i- u: e! r8 U    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

* j% E' Z' L" u  X出错 Untitled2 (第 20 行)
/ K" a% u# Z( @& ?5 Y+ M    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);


% H( v$ \) D' {3 G感谢大神们！
# }# c& k& |' V% E7 N9 t# W) k

llbnmo

clear all
clc
. ?! v# D* U- g3 G$ `warning off

eta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
G1=2;   % [Pa]
" \  {( n# l* T3 Q0 NG2=5;   %[Pa]
  L6 J6 ?7 i- a$ Qyieldstress=5;  % [Pa]
. h0 p$ Z+ L3 G6 S' Wcontactstrain=5;
; P; U, c# e; i& O

" N" [/ F) A, {$ c; Q: A) H7 r- _flow='creep & recovery';
( B: g$ ~' T2 O7 x8 f& `stress0=1:2:10;
$ F+ ~# L; h& D) C- \risetime=0.01;
$ F: x" h; z% K% J5 screeptime=100;
for i=1:length(stress0)
* ]( l) p! E3 x7 ^. ^$ v- z( `& }    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
0 N4 e* c! ~8 `: K! j    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcat(num2str(stress0(i)),'(Pa)');
end
xlim([1e-4 1000]);
xlabel('Time(s)');      ylabel('Shear strain(-)');
& a2 r9 \- S3 T$ e& n. flegend(a)
# C' ]4 U" e: y" Y: {: K) k+ Lfunction s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
    case 'creep & recovery'
. n& O: u  `6 l        stress0=arg(1);
5 m7 ^& i) `  Q& C: r! \/ O        risetime=arg(2);
        creeptime=arg(3);
        if t<risetime
            stressrate=stress0/risetime;
* V! e' l1 \9 L3 ~  j; A. N            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
2 a5 \1 L% i; o1 V            stress=stress0;
        elseif t<creeptime+risetime
& Y& `+ d" [" B+ R            stressrate=stress0/risetime;
5 X! _: l% O, F            stress=stress0-stressrate*(t-creeptime);
        else
; I1 ~  W' I1 i2 ]% S, u5 T            stressrate=0;
            stress=0;
        end
        gam1=y(1);
  N. C# u' o5 o& a. I  W+ N- Y. e6 p        gam2=y(2);
5 A" F9 Z% }6 o- G$ c        s(1)=gam2;
2 i3 o; ?6 c  }. q) {        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);
1 T7 h! t8 s5 A5 K& P- ~$ ^        s=s';
end
5 V, J+ N/ c  M' B* Tend

小白的白

楼上正解！

IBB-EUT

来学习一下