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

cichishia

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

% ^3 f( M, o& [8 }% a& b  E. ?这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
case 'creep & recovery'
( L/ ^0 B% O! o        stress0=arg(1);
5 J# a% f! Z" s3 L! t! m        risetime=arg(2);
% c8 n  z* A( h' I0 B2 W2 X! x5 W        creeptime=arg(3);
        if t<risetime
( A. A1 G: I$ H# B: [9 \            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
7 O+ j' g1 b7 l            stress=stress0;
( ~8 K7 a$ Q0 @7 e/ r3 a# F5 x        elseif t<creeptime+risetime
% V; N; J4 F2 i) z" J2 g            stressrate=stress0/risetime;
! c3 X& @7 S' d' f* s' x! V            stress=stress0-stressrate*(t-creeptime);
        else
5 b/ E+ G9 c( u& X            stressrate=0;
( @" F; e5 F* ~3 M  z            stress=0;
        end
        gam1=y(1);
& l. b* E! y" F- a        gam2=y(2);
        s(1)=gam2;
6 c% T  A4 l! P/ N- n3 h        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);

/ \; z- G, F: S2 T1 ~6 q# t

1 _3 D( p. P# ^, b8 Q' _clear all
clc
; e2 F) f( K+ X: i. Nwarning off
" `2 ?$ {+ H5 c+ d& ?. P/ Z9 y/ `
, z, i  ~+ `9 G3 ]- O2 xeta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
G1=2;   % [Pa]
G2=5;   %[Pa]
yieldstress=5;  % [Pa]
: ?2 r1 |5 X( P7 ?3 p$ W0 scontactstrain=5;
4 h, L$ y' H; r6 }5 n, J8 Z! _

flow='creep & recovery';
; m# s/ v) j& M& |# c5 x' dstress0=1:2:10;
5 q; ^: ]9 c7 l, K* J6 prisetime=0.01;
creeptime=100;
for i=1:length(stress0)
) {6 j* R1 x1 l    [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)');
end
4 I+ |) Z3 U/ t" vxlin([1e-4 1000]);   
( `5 _6 G5 f7 pxlabel('Time(s)');      ylabel('Shear strain(-)');
" Z& t( k' G' t7 u$ {' V7 f# Q' wlegend(a)
, w; E. F& x( S7 l' K
- l( }$ c+ u9 H7 r2 C
* Y' N2 m: R1 f: ?. I
3 x( S4 U" M! L错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。
# g3 u0 H; c# x/ i- p# s
出错 ode15s (第 152 行)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
. }5 b& E$ Q, ^$ h
出错 Untitled2 (第 20 行)
    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);


! a' L( ~' i- R感谢大神们！
  ^7 n5 z8 ^% r8 K! W- x0 \0 m9 W

llbnmo

clear all
* a. w/ y% R2 v( }: Kclc
3 t" o7 [  d- ~# j, }4 z% Uwarning off

& g# ?3 B+ A: b, c4 p* x! Neta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
1 M0 \& o; J- l4 v. [/ G9 AG1=2;   % [Pa]
9 s/ O" _! B' J6 g' qG2=5;   %[Pa]
* @' q7 Q6 O1 b5 e+ n# _2 \3 [yieldstress=5;  % [Pa]
# O7 ^9 R- R" v# r, m6 acontactstrain=5;

' [) f# ?' [7 n# l2 A7 `. K$ \
flow='creep & recovery';
- W( e3 `6 w+ n  Y+ X7 i) d1 Kstress0=1:2:10;
risetime=0.01;
) {% q# y! w) K% S$ `/ fcreeptime=100;
7 p5 C& u4 r% Hfor i=1:length(stress0)
( {% i% O7 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}=strcat(num2str(stress0(i)),'(Pa)');
+ y# C9 l7 ?5 m7 Y$ y2 hend
! v- b$ X. |7 O; q: i5 |xlim([1e-4 1000]);
  e* k; F7 i/ b2 v( Axlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
    case 'creep & recovery'
        stress0=arg(1);
        risetime=arg(2);
% {, V# m# r+ S1 {        creeptime=arg(3);
( L: A* }) w* u, d; I        if t<risetime
1 K) P: j8 b$ S' @) \+ W            stressrate=stress0/risetime;
2 \2 a* V; ?; e            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
" M- p% R# h0 p9 e# B            stress=stress0;
        elseif t<creeptime+risetime
            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
, |  E) _+ C* s4 G4 R) r4 E        else
, ]1 O! @0 O6 I1 L5 g            stressrate=0;
            stress=0;
0 w! z: c; Q/ Q, _: ~6 u! r        end
0 O# C5 \) |% S( X9 i# O" `) ^5 v        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);
        s=s';
end
end

小白的白

楼上正解！

IBB-EUT

来学习一下