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

cichishia

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
: p' E: l% @; X3 J7 ^8 U
这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
$ c, W' Z; L7 L$ wswitch flow
6 @1 |  U& ~& j. s4 x6 _case 'creep & recovery'
' M* O' y8 @1 `& X/ g3 k; {        stress0=arg(1);
9 g- d6 i" d+ f1 ]; ^( z6 Q        risetime=arg(2);
        creeptime=arg(3);
) D7 J5 y: a; f$ k. P! O* F        if t<risetime
            stressrate=stress0/risetime;
9 x. H. u1 f: j            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
5 B/ p2 e2 E+ s0 f6 b& D+ w            stress=stress0;
        elseif t<creeptime+risetime
6 V/ r. \$ Y3 ~; p8 C) w3 f: {            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
$ ]8 g# Q8 z# I# o- s        else
            stressrate=0;
            stress=0;
        end
        gam1=y(1);
        gam2=y(2);
. u* _" Y! b! E3 E1 u5 Y        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);
- i! p0 _9 a; T, w7 e" g

, e; m- |! z  x
* \; S5 M" A1 z* e/ `clear all
: [  P' k2 s) o: o# Gclc
+ j3 w" ]! D/ s' W9 g$ L1 [- z5 Q# kwarning off
2 n, L* w) n. ?$ B8 @/ O( o
, p; \: s1 y- C4 ~* h  x0 J* c* y' r4 Ueta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
$ U9 G* f7 {, p0 W" l5 g/ E2 B6 uG1=2;   % [Pa]
G2=5;   %[Pa]
* o4 ]& p# s4 ^; e0 byieldstress=5;  % [Pa]
contactstrain=5;
5 f3 P8 d9 G$ \2 l! p+ g& s

flow='creep & recovery';
, ]% h% i) C  Kstress0=1:2:10;
risetime=0.01;
' G8 U  @/ q  h# J2 X$ B- @+ Ucreeptime=100;
for 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
    a{i}=strcatnum2str(stress0(i),'(Pa)');
" @/ S; s7 V, e& Cend
xlin([1e-4 1000]);   
- W3 ?# q6 H1 F; b- zxlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
* ^; T$ j# ?. c7 O6 K* Y


# N# i! u" j1 a9 X) ]3 `* {& p错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。

/ h, a/ T6 u5 I$ m- Q出错 ode15s (第 152 行)
# L' Y; |/ ]5 x' T    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
: J: {) R& x8 J' `( D) j& Z1 L% m
出错 Untitled2 (第 20 行)
    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
9 {% i4 t- m/ p2 M
* l$ y" M( Z. L, Q, e+ G! k
3 `0 A# J7 {2 f, |感谢大神们！

llbnmo

clear all
clc
warning off
$ a# R& ?# O5 D& i
# k' f0 g" o5 [! {2 I& ?( ueta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
" f+ S& H$ F1 R% J4 k6 @6 `* ~$ ?' {0 U, uG1=2;   % [Pa]
G2=5;   %[Pa]
9 Y+ \; y& `6 X5 ^$ jyieldstress=5;  % [Pa]
6 v/ {) y! n: }4 Gcontactstrain=5;

/ z6 y" [* _- O5 f% E0 W2 }& Y
flow='creep & recovery';
, S4 T2 H: E+ X  S  rstress0=1:2:10;
risetime=0.01;
8 s* x1 |, I2 m& |( ]0 |  _& hcreeptime=100;
for i=1:length(stress0)
  u1 I8 e, `) ?, g/ H    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
1 x9 \# R1 m7 _8 k" l' ~( J    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
/ U. u' C# j* V0 W$ |, O& J0 ^    a{i}=strcat(num2str(stress0(i)),'(Pa)');
end
6 @3 f5 s- x9 J3 Z" M$ Kxlim([1e-4 1000]);
xlabel('Time(s)');      ylabel('Shear strain(-)');
5 x) D! G/ M0 {1 G! Z& g( Vlegend(a)
, \1 D& z: C- bfunction s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
5 z& T! W/ w/ c: J3 s2 B2 T! bswitch flow
1 R+ }1 W6 @7 _7 U    case 'creep & recovery'
        stress0=arg(1);
  @! O  d2 I4 H: F        risetime=arg(2);
        creeptime=arg(3);
/ s1 r4 X7 y2 [        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
  @, a9 G+ P9 O6 b6 ]/ L        elseif t<creeptime
            stressrate=0;
            stress=stress0;
& m7 r; B+ i: z9 b0 L        elseif t<creeptime+risetime
# X; h5 Y7 A( Q- G; E            stressrate=stress0/risetime;
) G) E3 _5 S; e! A& B" |. a7 L            stress=stress0-stressrate*(t-creeptime);
        else
" k4 X* v# Z  s9 I  O            stressrate=0;
  b6 ]# H4 W% ]" a, u            stress=0;
' S6 Z( k* c4 D        end
1 V; d2 b4 R2 V8 _) n! ]$ ]        gam1=y(1);
' j/ o8 w5 k' g        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);
; i( t" i/ I( u% W4 s3 F. {/ i( j& t        s=s';
2 Y' ], ^9 ?" w' Bend
8 X# t1 H' E  C0 vend

小白的白

楼上正解！

IBB-EUT

来学习一下