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

cichishia

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

6 y2 Y9 A- _5 O% \3 k# H这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
  P& Q+ W5 N/ I9 d! Z8 |! o% Ncase 'creep & recovery'
) ^0 Y7 ~& ?1 B% U# O+ h1 P        stress0=arg(1);
        risetime=arg(2);
/ m0 o( V  R/ l: k- f3 v3 Z8 {0 g1 V! A        creeptime=arg(3);
        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
, }6 ~! X) n/ U8 ?  C. Y            stress=stress0;
        elseif t<creeptime+risetime
2 ?& a' ^9 q% l- j1 B3 R- s            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
4 I; f' F! Q( b' y9 B4 V% W            stressrate=0;
            stress=0;
! K# i) W7 ?2 U+ C        end
0 E4 O2 O7 ~7 A& p2 c- |        gam1=y(1);
8 M4 h8 e1 `% ?3 e+ P6 A: E6 R        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);

( ~5 y% z7 A7 K  n# m# Z
4 P; L9 R8 |; K2 ]2 S$ U; D: Q9 |6 t" o. a
clear all
# d; _6 \5 ]; |( o% u9 U6 h$ X; _clc
warning off

1 Z- D$ u) y5 ^" ceta1=10;  % [Pa.s]
! B7 M7 D0 e- v9 n# H& }1 oeta2=20;  % Pa.s]
G1=2;   % [Pa]
G2=5;   %[Pa]
yieldstress=5;  % [Pa]
contactstrain=5;
3 z4 r# s2 m: t/ X% @
1 b$ B' u8 w6 P- A9 Q7 ^* q
flow='creep & recovery';
& F' n, _% [6 {; fstress0=1:2:10;
# h# H6 _! i2 o1 r7 w& ?& h" jrisetime=0.01;
creeptime=100;
for i=1:length(stress0)
" m/ I+ E8 ^4 G/ g. O6 I* L9 {    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
& v8 o1 g2 \% m0 F8 a( f( E5 u$ ?5 U    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
1 ?) _4 y- P* Z, Y; {/ e    a{i}=strcatnum2str(stress0(i),'(Pa)');
4 F) E0 T& _: [+ G1 F! _8 B6 U! @* xend
; H! r; j  ^$ I+ d; _xlin([1e-4 1000]);   
xlabel('Time(s)');      ylabel('Shear strain(-)');
9 }; K( c( f$ Z4 @5 }8 _legend(a)
( f- C1 U4 i' q7 j) ^* B

, q* n0 G+ t" N7 J* G3 B) ~
错误使用 odearguments (第 93 行)
, |1 l) K. g. ~# D; `9 [( bMYWORK 必须返回列向量。
, `0 w4 N4 J7 _4 {* X- O+ t3 o
出错 ode15s (第 152 行)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

$ k5 H$ X' U2 ~! t* A; p6 V出错 Untitled2 (第 20 行)
    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);


9 I" Y4 q7 ^) s: }感谢大神们！
* P* W, W- D1 c4 Q; y! l8 [# T

llbnmo

clear all
' V( j. J$ y1 G' W3 `- }+ \clc
' R6 q8 U) ]/ u1 S; fwarning off

eta1=10;  % [Pa.s]
3 F7 u# s- H# aeta2=20;  % Pa.s]
6 T2 X% O) c. A5 L& d. aG1=2;   % [Pa]
5 s! c4 p+ N3 X3 Y% n, IG2=5;   %[Pa]
yieldstress=5;  % [Pa]
contactstrain=5;
; n1 v- o" p: R% {1 S

flow='creep & recovery';
stress0=1:2:10;
" w) W6 p+ e; u; U1 G6 erisetime=0.01;
" P$ n' Q# w  p4 C* mcreeptime=100;
for i=1:length(stress0)
    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
' u3 ]. r2 C$ a! [' g9 q* f    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
' R3 h6 ]) i6 l0 V0 m& A  j    a{i}=strcat(num2str(stress0(i)),'(Pa)');
end
( Z. _* q) f& Mxlim([1e-4 1000]);
. s/ p6 h6 G2 F- d( E" dxlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
+ F1 C( [) q- |, g! A) m& F8 M5 S    case 'creep & recovery'
/ A! S% H% M* I        stress0=arg(1);
        risetime=arg(2);
        creeptime=arg(3);
; ]7 ^: T1 I/ U9 @$ d' x( @        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
1 V4 y: T) o! I, @  n+ u' R            stress=stress0;
        elseif t<creeptime+risetime
            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
$ u6 c+ F$ g9 W# R8 C            stressrate=0;
8 \- E( G7 N4 D1 R            stress=0;
. l6 B  p& ~6 n2 D" d7 b$ o) r, f9 u        end
        gam1=y(1);
1 x( ]  \% p# Z) X% C% b        gam2=y(2);
& _) Z- Z3 m- ^* W        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
0 Q  B# V- f/ T. f  o: w/ @end

小白的白

楼上正解！

IBB-EUT

来学习一下