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

cichishia

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
/ y/ g8 Q+ \& s' ~! R
6 d2 z$ F2 k# _  t2 a/ k) s这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
case 'creep & recovery'
, k' W% ^! {0 m: B: d  s# S( ?        stress0=arg(1);
        risetime=arg(2);
* @/ `- F' x! p+ A; d/ Y3 e% X        creeptime=arg(3);
        if t<risetime
- Q% D  N$ U& ]$ q9 P            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
0 D' A  V8 R2 z            stressrate=0;
" H4 Z2 Y) U  I$ U            stress=stress0;
        elseif t<creeptime+risetime
3 b5 p2 L$ O3 ?9 `! [4 l+ J# z            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
: w2 g' M/ i% T& t- O6 ^$ E# @            stressrate=0;
            stress=0;
        end
        gam1=y(1);
& D. J- b. M$ ?+ l2 \& b9 ]- p        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);
1 r! w+ p& d% p$ |
% U, o9 ~. _5 T( t# j6 D; ^2 d  u" _

2 E+ g9 h$ x4 D; _$ oclear all
( z4 E! `# L$ C2 T% [1 Pclc
warning off
; y9 X9 H0 i1 z" @2 M) O- k
/ v/ e. Y# P1 O( |* ~8 Beta1=10;  % [Pa.s]
+ X/ {. s1 k+ r) P4 d$ E' K3 feta2=20;  % Pa.s]
- i% W" c6 `+ k- ^1 Q  m5 SG1=2;   % [Pa]
8 B2 s& |9 H: D" T( K+ nG2=5;   %[Pa]
yieldstress=5;  % [Pa]
; n9 b: m5 K% I$ J' @* Zcontactstrain=5;


4 v$ w# t% U$ E6 t9 Jflow='creep & recovery';
stress0=1:2:10;
6 v; C9 C+ C+ xrisetime=0.01;
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]);
$ D& w, \! `5 ^- o+ m: I' l9 s! G    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcatnum2str(stress0(i),'(Pa)');
* C$ L) L1 H+ K+ eend
xlin([1e-4 1000]);   
xlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)


: ^) l6 p0 E0 o$ }
& l9 U8 Q2 t0 y/ u# r错误使用 odearguments (第 93 行)
4 v+ V4 I3 F8 u: B5 \MYWORK 必须返回列向量。

3 r& b' o. H3 V/ g9 l- v6 Y出错 ode15s (第 152 行)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
( p$ _1 {/ {2 T1 _) U8 s: h* Z: O) ^
出错 Untitled2 (第 20 行)
    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
; {9 e5 y3 K6 }  w. a2 Y

; J# v* ?6 h+ @! d6 l; Y感谢大神们！
5 Y: l8 Z1 w) |( W+ D1 C

IBB-EUT

来学习一下

小白的白

楼上正解！

llbnmo

clear all
clc
warning off
) [" v& M( Y, v' F5 C' L. v
9 `7 Q9 E* a3 @8 `  y# reta1=10;  % [Pa.s]
* X! E+ Y& w- n6 H3 F. i: Keta2=20;  % Pa.s]
G1=2;   % [Pa]
G2=5;   %[Pa]
  F  Y- f3 m  B0 D; S3 _yieldstress=5;  % [Pa]
' c6 C5 `# T9 J$ H) c7 W5 m5 U; X3 Econtactstrain=5;

8 w4 G6 ?+ m) ^/ C8 y
/ ?+ D! n9 ~- M( c! Rflow='creep & recovery';
% j5 w) l  {7 Y0 F% Sstress0=1:2:10;
. S8 D- I8 a, G* H1 V& r' [risetime=0.01;
. d$ j; D- v/ w* [9 {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]);
" F. `8 Y% Y1 v1 ?    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
    a{i}=strcat(num2str(stress0(i)),'(Pa)');
' {; W6 a% u. C5 C& i8 J' dend
xlim([1e-4 1000]);
xlabel('Time(s)');      ylabel('Shear strain(-)');
legend(a)
( Y! e- H' u6 A# Kfunction s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
: _7 w  S) G/ e. {) q/ Lswitch flow
    case 'creep & recovery'
        stress0=arg(1);
/ F/ T% I% t1 s; f* X9 e" N" R0 `        risetime=arg(2);
0 |0 G2 T( t6 _2 L, V. c0 ?/ v! Y5 I. z        creeptime=arg(3);
" T* Y' x* u# ~$ o        if t<risetime
1 h( }# g- I+ i            stressrate=stress0/risetime;
            stress=stressrate*t;
        elseif t<creeptime
) h" m. t! N2 }            stressrate=0;
& T$ h6 Q9 L. A+ l7 B: _- b0 h            stress=stress0;
        elseif t<creeptime+risetime
3 y% E9 i# h; @: X. z7 y            stressrate=stress0/risetime;
            stress=stress0-stressrate*(t-creeptime);
        else
            stressrate=0;
            stress=0;
        end
2 r1 L( b5 Q/ f: k1 e1 r- m3 X" q: G        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$ Y7 J7 }. J8 m; f3 d        s=s';
end
end