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

cichishia

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
  ^" j/ f% a) i" C) R1 b$ n: m1 [
这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
: P! W7 [4 E6 ^6 M1 h- N$ @' ]# v, Bswitch flow
( ?! z$ J' `2 T# }7 }case 'creep & recovery'
        stress0=arg(1);
1 {) M" w& Q; U( v* F9 k' b6 N" g        risetime=arg(2);
+ k# ?) S  V- z/ p4 P  m        creeptime=arg(3);
0 T4 S) M' \; Z/ f5 ~8 c- _        if t<risetime
            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
            stressrate=0;
            stress=stress0;
        elseif t<creeptime+risetime
" Z6 K( P/ E' S, C. L" b3 p            stressrate=stress0/risetime;
: ~8 L1 H4 i7 m            stress=stress0-stressrate*(t-creeptime);
( T& n; U0 W" c% J; N- r# q' d        else
9 o  v$ }! M3 g7 p$ Q$ L% N, H            stressrate=0;
            stress=0;
6 \6 y0 Q1 c" B; C        end
        gam1=y(1);
, J6 s4 t* l/ ~3 I- V! X- H        gam2=y(2);
        s(1)=gam2;
+ \. d5 i3 o$ P3 P& J9 ^        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);


; U5 R) v7 J/ {8 A
clear all
/ p: G: f4 P% ~5 i* jclc
/ O. n" t  ~7 L$ s0 w% Awarning off
' v0 f8 m* U( h/ ^  `# C8 @
& Z2 X' O" {1 o' X: U* Beta1=10;  % [Pa.s]
- P, g' w: D1 T: v0 h. _eta2=20;  % Pa.s]
G1=2;   % [Pa]
G2=5;   %[Pa]
yieldstress=5;  % [Pa]
1 o% u  j4 |/ }contactstrain=5;
" o, ]3 I1 W% T. v7 W' p& u
1 d: ^& N/ `" U0 E
flow='creep & recovery';
/ r  X4 j  g- S$ f0 ~+ estress0=1:2:10;
: z0 C1 R3 h3 l/ xrisetime=0.01;
% M3 o$ F( ~: \0 M+ H4 ocreeptime=100;
for i=1:length(stress0)
" ~; B& F3 S% x+ n7 R9 X- b1 V    [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
7 y3 L$ i- H1 S( m" D    a{i}=strcatnum2str(stress0(i),'(Pa)');
) t- D- q1 e0 X5 o. Iend
xlin([1e-4 1000]);   
( n6 S; o/ I, ]9 n0 n; M! Zxlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)

  b. T6 b  J$ }

错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。
/ K  T* \$ b. _
出错 ode15s (第 152 行)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

0 o5 W6 j: T3 ?  X. N/ g; I出错 Untitled2 (第 20 行)
# E$ X! I* |, H  l5 T1 C- @1 }    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);


( P2 ?3 `% M4 M: W感谢大神们！
2 @" ?$ v3 ^; O& N2 u

llbnmo

clear all
" x" [1 e' G8 I9 S$ u! L2 }, bclc
/ X/ L- a5 l2 U- x& |' fwarning off

/ W6 D& n5 [7 E- Teta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
3 r8 a  r/ `1 P  w! zG1=2;   % [Pa]
* Q4 U/ @# N  B6 MG2=5;   %[Pa]
" E. ?* S! q* w+ p* h+ ryieldstress=5;  % [Pa]
contactstrain=5;
" R$ e; F9 j( d: R5 B- }! L# k

flow='creep & recovery';
stress0=1:2:10;
risetime=0.01;
9 V: w9 w! O0 W- n0 p4 @creeptime=100;
for i=1:length(stress0)
9 \. x! V; g2 e9 ~9 k, H    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
1 e" F; J8 }/ l1 V# [+ M; l% c% q    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcat(num2str(stress0(i)),'(Pa)');
$ w6 D+ D* t% j  u# T: nend
xlim([1e-4 1000]);
xlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
7 G: o* Y8 O  s% M  m, Kswitch flow
    case 'creep & recovery'
        stress0=arg(1);
  F% D6 H9 b# |- ^5 W2 d0 L        risetime=arg(2);
/ x+ b+ P1 z; r. ~/ ]        creeptime=arg(3);
        if t<risetime
0 ~1 T) {/ d! ~& Z; d  I& m            stressrate=stress0/risetime;
" B' F  {% W+ ~; Q' i& x            stress=stressrate*t;
1 q8 G% R9 Y8 J        elseif t<creeptime
4 @3 ~& J& z: L4 M" w7 {            stressrate=0;
) w. h- |* P" B0 Z- t            stress=stress0;
# v# h) G( O/ w& p) l        elseif t<creeptime+risetime
; u0 L  Z' |- E* i  L            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
            stressrate=0;
8 l* H: W4 I6 T0 j) Q            stress=0;
+ N6 O/ `0 y( x        end
        gam1=y(1);
        gam2=y(2);
        s(1)=gam2;
: ~7 i7 t+ @0 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);
        s=s';
: Y6 G! t% ~) \5 s8 I* Gend
: q2 F% C9 V7 Y! m7 X* G1 e+ nend

小白的白

楼上正解！

IBB-EUT

来学习一下