EDA365电子论坛网
标题:
matlab解常微分一直报错返回列向量
[打印本页]
作者:
cichishia
时间:
2021-1-5 10:59
标题:
matlab解常微分一直报错返回列向量
我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
0 x7 ~+ A" @5 C: ~
3 E! g3 u' v& `9 u. N& Y9 o8 Z
这是代码：
~5 M# _3 Y0 m! K6 t7 Y
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
: \5 M5 _8 j* x2 i) V* Q
switch flow
# i. s2 k0 \$ I' l9 I9 @
case 'creep & recovery'
- i X" f$ e3 Z+ m. w4 h5 _
stress0=arg(1);
( P& G; A4 K3 i1 x( G/ l
risetime=arg(2);
. o7 S7 p" h% i( H
creeptime=arg(3);
2 G( L: P1 Q7 e- }% U& j w7 W# r$ {
if t<risetime
+ i" }) p- M- d$ m0 n) u
stressrate=stress0/risetime;
7 u. {5 Z7 a% G8 K$ N3 C6 K0 C
stress=stressrate*t;
* C% k' L* U; d5 G, B @
elseif t<creeptime
: _/ s; y. _0 ]: \: f* E* n) h
stressrate=0;
- U9 D/ }* s/ p- T; @8 E; N
stress=stress0;
0 f8 e+ a0 i* k1 Y) M7 Z
elseif t<creeptime+risetime
7 T* {& z$ @/ q' s1 a) i
stressrate=stress0/risetime;
" ?6 ?" ?" c8 `1 Y4 Y5 D. ?- j
stress=stress0-stressrate*(t-creeptime);
+ d: B. w' T' o$ D9 p9 }
else
7 T* w1 h; R6 y( ~* U
stressrate=0;
7 S# {0 X4 b* U8 W# O% b9 [* G
stress=0;
; d- v1 r7 }+ Y. q# ]% F: r
end
- O0 f$ k) a! B8 X
gam1=y(1);
! D5 j. c) E7 B) S: y, k
gam2=y(2);
5 a. p5 j' n9 b; s0 k
s(1)=gam2;
" W% j) ?/ `/ r7 Y+ }9 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);
( M% K9 ^' Y$ P0 \
$ S0 M$ U& t' F; C
8 I1 [6 M) c7 u8 s
( m8 A3 j2 i0 ^- Z, f$ M, n
clear all
% x5 G) l% q8 M4 L" h l
clc
3 z: ~4 `4 Y6 L* o4 U, [9 A) |5 Y
warning off
' q" L' `4 b- H& L0 K' Q! c
8 e! y6 H9 I2 x. W1 e! K0 ^% h
eta1=10; % [Pa.s]
" H# c$ V) e4 h6 }% Y1 w' E
eta2=20; % Pa.s]
+ f2 ~: y7 s0 b* w! h* v! W
G1=2; % [Pa]
5 h3 y- z3 ?# ^4 e
G2=5; %[Pa]
' s* v8 n$ k- Y3 c* h5 A
yieldstress=5; % [Pa]
@. f' C* A. Z0 U8 H
contactstrain=5;
2 a8 W. Y3 f2 F
) X6 ^" n" m) P8 k7 W2 k& o3 G
( W) l# \4 Q1 u! D- J4 D! J$ o# g
flow='creep & recovery';
% T9 P& Z7 ], e+ j
stress0=1:2:10;
; D# l2 F. J/ x" b% `) |& H% r
risetime=0.01;
" r- g2 `+ O0 j* `1 A y. i
creeptime=100;
' O2 g! ?( h2 A+ F
for i=1:length(stress0)
5 Y2 s$ \; H: X, |4 B0 p
[t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
% V/ s0 K% m$ x* f! \" U' L1 c
subplot(2,3,3); semilogx(t,strain(:,1)); hold on
* s3 B2 Y3 P4 a) \5 M
a{i}=strcatnum2str(stress0(i),'(Pa)');
. ^ O7 B9 q; T# T3 f
end
( K. R; i8 W( c0 j5 W. s/ H3 F- H
xlin([1e-4 1000]);
7 \3 W. y& _% H- V& F% I
xlabel('Time(s)'); ylabel('Shear strain(-)');
1 U" m1 A( \' T/ {- X
legend(a)
! f# U4 j- J, r. t2 G+ p* f
: P4 p# ]# F5 B; v; S6 b
6 B. O4 T) a: ?1 z* a( b! L# M, z
, Q0 f7 _. \/ K9 {* B9 W, ~3 t+ ?. ^
错误使用 odearguments (第 93 行)
, C8 |* ], _5 T
MYWORK 必须返回列向量。
1 W$ h* @8 H4 ^1 t! z! x
& _& K+ i/ v5 j: g0 I
出错 ode15s (第 152 行)
4 f- I& |* _* g1 H- }
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
9 }8 u r3 H# [5 C V" t
\4 Z* l/ t; v$ |# o1 [
出错 Untitled2 (第 20 行)
! P$ F3 D! x$ u/ {$ Y
[t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
% G: T2 e4 @' V9 O r5 j
: t* F( y1 W- A5 ^2 W `$ b* J
( w' Z* o; C$ q7 ~. r& o2 W
感谢大神们！
& M) E5 ~5 d4 l* u. _) ~
作者:
llbnmo
时间:
2021-1-5 13:09
clear all
' ~9 L1 _ ?3 f; Y
clc
1 ?. w1 ^3 q3 K0 l+ X
warning off
) N' Z7 D. o1 k/ q$ m8 V0 n
' |- ~! I* G% J# l2 ?
eta1=10; % [Pa.s]
% c2 b7 v# R3 q
eta2=20; % Pa.s]
/ i T* W4 f% d
G1=2; % [Pa]
( n' [4 _5 H6 h. g
G2=5; %[Pa]
! c( R q2 |" o0 o, @
yieldstress=5; % [Pa]
; j' ^: h' h- w w
contactstrain=5;
& |4 {4 }! z6 U3 g& L+ ]/ W( i
6 z9 \2 N+ n/ i) w0 D' O
- K: T5 Q( L9 f, Y
flow='creep & recovery';
+ ~- B: y2 x/ l# h* r
stress0=1:2:10;
5 n! ?& \1 ^7 l. p% Q$ g
risetime=0.01;
3 i2 t5 \+ R$ G* i+ v4 c
creeptime=100;
6 G0 f' W& ]; ] D' W
for i=1:length(stress0)
# M9 k$ G9 m3 E
[t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
0 J$ |6 f$ i" G, n4 y
subplot(2,3,3); semilogx(t,strain(:,1)); hold on
8 ]4 T% w9 j( X
a{i}=strcat(num2str(stress0(i)),'(Pa)');
6 K8 ^7 v, z5 K
end
. h3 B( d4 b; ?4 B- K" O6 Q
xlim([1e-4 1000]);
6 ?, R# j8 P& `! f
xlabel('Time(s)'); ylabel('Shear strain(-)');
+ ?7 e5 V3 v0 }5 z5 ^
legend(a)
3 ?% j5 i4 A- y
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
3 [6 c& K0 ~0 u
switch flow
7 E0 ]1 _* |2 r
case 'creep & recovery'
/ R. f) m; B( n
stress0=arg(1);
# \! n4 I% K& H" X- b6 X; K
risetime=arg(2);
7 U7 S5 y& T9 @/ Y. x
creeptime=arg(3);
. r" y. X" `/ a
if t<risetime
3 K5 R4 z9 K& k1 J. F s- [
stressrate=stress0/risetime;
- _4 Y; l" {# Y+ O) O" P
stress=stressrate*t;
7 R* g& F* @& P( F/ r7 x. i
elseif t<creeptime
' G( w6 _- D, `/ @
stressrate=0;
6 V: D; _0 H# E- i
stress=stress0;
8 v) E3 i. i& j
elseif t<creeptime+risetime
9 }( O# i! g" I# Z6 ?
stressrate=stress0/risetime;
/ k, g# ~; K7 I
stress=stress0-stressrate*(t-creeptime);
( A& _* [9 e* `
else
/ a0 v' @+ v* ?* }! c
stressrate=0;
. U) w3 U5 K5 J3 a
stress=0;
8 `' [" m; \& p9 C
end
- m l7 p3 w8 Q
gam1=y(1);
( [- O3 \+ R" b0 x
gam2=y(2);
9 a2 V3 z' P; q, ]% | E
s(1)=gam2;
/ e' j& V& y' W- N: G- ^
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);
5 I- p, X% Q: j5 L! j( R k- h/ W5 B
s=s';
( @! Q9 z2 E1 {! q5 g% z
end
% y( B, i$ Y G
end
作者:
小白的白
时间:
2021-1-5 13:17
楼上正解！
作者:
IBB-EUT
时间:
2021-1-5 15:43
来学习一下
欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)
Powered by Discuz! X3.2