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

cichishia

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
( S2 F+ ]7 q4 [  r0 |
这是代码：
- I8 s$ i9 y# u2 m3 Dfunction s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
case 'creep & recovery'
0 e  y# l; D6 R8 q. a0 Z. T( N% P1 v        stress0=arg(1);
9 b! \: A3 s. y2 x3 U        risetime=arg(2);
# n" h. f9 h6 B7 s( N; g% M& u        creeptime=arg(3);
( N9 O1 B( k& `4 b5 j# H- m        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
, n# _# T3 f$ \8 n            stressrate=0;
            stress=stress0;
# ]% s6 T2 e- N% y        elseif t<creeptime+risetime
            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
            stressrate=0;
$ R% I! y4 k8 g$ B+ b. B4 i            stress=0;
        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);

( T  O) {2 v3 N3 O5 Q3 q+ E

5 p% M; l/ k" X9 r) `clear all
clc
warning off

; j7 Q/ b: ^, x0 v5 |eta1=10;  % [Pa.s]
/ E( Z: z# y0 |! ^3 [* U4 q; ueta2=20;  % Pa.s]
3 \! G. g. s9 uG1=2;   % [Pa]
G2=5;   %[Pa]
$ d# o% ~$ g; E7 v6 ]* Zyieldstress=5;  % [Pa]
  z+ U1 o" }: F7 w. Z* ^2 Fcontactstrain=5;
4 n( e/ U6 z; A2 t2 I* T" T, C

flow='creep & recovery';
1 Z5 w& S# z; `stress0=1:2:10;
risetime=0.01;
creeptime=100;
7 [* z' P3 ]( h, h/ B' p# Gfor i=1:length(stress0)
    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
: G9 p( x9 Z7 q9 b" G# H) T    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcatnum2str(stress0(i),'(Pa)');
( a1 {6 R9 b; V# d4 n$ A! W9 q0 Tend
/ A! Y$ [% U" }xlin([1e-4 1000]);   
xlabel('Time(s)');      ylabel('Shear strain(-)');
6 M# O% s9 T! i$ Flegend(a)
, k. r- {( d( l' a! d
+ \% E4 Q8 ^% q7 M
# @: I$ m$ V7 W, A* }$ G6 t
错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。

$ t/ ]% k$ ^6 f( V. }( w2 w出错 ode15s (第 152 行)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

( W3 |# l8 {* b) O& o& c- N出错 Untitled2 (第 20 行)
    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
8 f1 ~: p$ s: K& E- [! j. a1 o
) D, O) K' F# o  P" s
感谢大神们！
, g9 y0 W0 i7 _- |" D

llbnmo

clear all
clc
warning off

2 B" `2 H2 }" c+ ^eta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
3 u+ ]8 l3 @6 `. c- x. ~3 D( u1 j  D0 DG1=2;   % [Pa]
G2=5;   %[Pa]
yieldstress=5;  % [Pa]
contactstrain=5;
" W" J% E; [2 T
' p+ ~0 R: g2 `; k9 [
) ]: C- \$ W0 r8 h$ k4 vflow='creep & recovery';
stress0=1:2:10;
) l3 V1 h- t$ o2 b& g; d0 C4 c% x& Qrisetime=0.01;
& d6 N5 g0 E# e9 K5 r5 \creeptime=100;
for i=1:length(stress0)
    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
8 y' c9 Z1 o1 s( ~9 j& @    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcat(num2str(stress0(i)),'(Pa)');
end
7 l) b1 d, e/ D+ k* k* G6 bxlim([1e-4 1000]);
) D4 G: w6 v, x, r. {$ hxlabel('Time(s)');      ylabel('Shear strain(-)');
+ w$ x1 }6 [" T& F- u0 O5 olegend(a)
2 a$ c% p) `# efunction s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
' r1 |) D, A$ g; X1 M, fswitch flow
    case 'creep & recovery'
" I: d9 {% k3 C5 I" X        stress0=arg(1);
+ E/ W( g- g/ l% J4 e        risetime=arg(2);
        creeptime=arg(3);
        if t<risetime
. S6 V) w. B0 }) D/ W- e            stressrate=stress0/risetime;
) g5 _% v+ F/ O/ \- l+ d            stress=stressrate*t;
) L3 V4 p. S7 V5 a4 }        elseif t<creeptime
            stressrate=0;
            stress=stress0;
        elseif t<creeptime+risetime
4 b; B, D: U, g1 e# A            stressrate=stress0/risetime;
' P* X0 ?+ j& s. A' A( r7 R6 f% }5 h            stress=stress0-stressrate*(t-creeptime);
- i  N9 n0 R5 d6 I; Y        else
/ i; m4 [! Y% x; `            stressrate=0;
+ D4 o  Z: K2 E3 h" L2 B- F            stress=0;
        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);
  H9 C# N' k! x1 a$ i        s=s';
- o# a9 S9 l: l3 @' Vend
end

小白的白

楼上正解！

IBB-EUT

来学习一下