|
|
EDA365欢迎您登录！
您需要 登录 才可以下载或查看，没有帐号？注册
x
我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛% R' [- D% l% C" w3 z4 }5 d
; H# t8 \8 X5 l1 J3 @2 a. g9 ?3 a
这是代码：; I+ l% }8 Q( i% A) s7 X* v
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
9 |$ e* Y9 T' ]5 q" uswitch flow- o9 }' Q9 `1 C( u* S
case 'creep & recovery'3 ` X7 h: L6 G1 ^3 A3 h4 P
stress0=arg(1);6 B6 ~5 Z, J X9 C4 j( H$ F) V
risetime=arg(2); v- [% b7 o: O) I* ?2 o7 `
creeptime=arg(3);
+ N1 G4 |; O2 }+ M+ R; ^+ _/ i if t<risetime2 c* a$ D- s" `8 i
stressrate=stress0/risetime;
7 A5 A u y( N& F. Z stress=stressrate*t;
, l0 V: L4 x6 K7 G elseif t<creeptime5 H; r7 V! }) s% l/ o" P$ h
stressrate=0;: Y+ S9 k" Y' `
stress=stress0;, O% L2 y# T; Y2 q6 ^' M
elseif t<creeptime+risetime
/ Y" p5 d0 x7 a0 V stressrate=stress0/risetime;
9 ?! ?; ]7 Y3 V0 c3 n9 H+ z; G4 r$ t4 | stress=stress0-stressrate*(t-creeptime);
! _. \: [7 n! Q( U else8 ~! g7 T8 K. a8 x$ D9 u* N o* n& S
stressrate=0;3 a5 l* w- R3 t, a1 V
stress=0;1 D( t, W1 S% ?+ A
end
* Q9 g& ^+ f# y gam1=y(1);
. D% s' j. }: `9 d! d4 X gam2=y(2);& s- J5 t3 Y7 P
s(1)=gam2;4 _* w: Z# R8 y8 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);
2 f( u5 h6 _# j7 a8 d5 r& c, l
& S* w: V: [0 z- f" w
( `0 k; H) r6 e# s/ Y7 u5 W% h. w# N2 h" M1 G+ T) e, o8 h" A% |
clear all" j, x& D8 C/ a+ @+ N$ \5 ^, t
clc7 a: V$ K, i7 X$ E' G( I( a
warning off6 |, k5 S# D. B. v
C% R; S H9 F: u& {+ Y8 m
eta1=10; % [Pa.s]
, r5 f7 [; L$ s oeta2=20; % Pa.s]# K2 j6 D' e( }8 y9 j; p) _
G1=2; % [Pa]! d4 y% h+ c8 l2 h6 \. @. d
G2=5; %[Pa]
, j0 y! v6 T: R# _5 [$ S/ ^yieldstress=5; % [Pa]+ C; Z8 O2 Z! j! s6 ~) F
contactstrain=5;
1 _# W- J4 \% F0 g% u+ F* K: j& b* h4 @) w \( g& g/ [
, n6 e2 i9 h! F8 ]flow='creep & recovery';! `- Q* C% t' ]0 v% c8 L/ ^" ?9 k6 w
stress0=1:2:10;
9 d% A' W4 w( f% T1 Q' `risetime=0.01;
1 e+ o+ Q+ F6 T- n9 ]) pcreeptime=100;1 W3 A% Q7 Y0 \+ b9 o* l
for i=1:length(stress0)
$ t' T+ V; h/ H! k [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
0 G+ z/ q; ~- k; T3 b- C$ u subplot(2,3,3); semilogx(t,strain(:,1)); hold on4 n2 z, v$ Y+ R2 L2 y4 `: ~' I
a{i}=strcatnum2str(stress0(i),'(Pa)');
- G" A/ e" k; e* ?# X& K7 Jend
; O9 }' N) `. _0 ^! Jxlin([1e-4 1000]);
* \& L+ Q/ t/ ixlabel('Time(s)'); ylabel('Shear strain(-)');" A8 _- n+ l6 j3 ?
legend(a)
! U* K/ n% }3 j
& Z( w: z/ E) u, s; T% U8 i# M0 t8 r# Z Z* W1 }; i* g# r
. ^" W' d1 H+ B0 Z错误使用 odearguments (第 93 行)5 R% y6 n" s% X/ |. a/ q
MYWORK 必须返回列向量。
0 w7 k- h4 C$ u. T1 k4 [# c
) @; A$ N: u% l+ Z. V- ^( {出错 ode15s (第 152 行)% g& @& c1 ~% S3 A6 [: U$ {
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
: H1 P7 Q1 E9 R1 Q9 M- D
- `7 t6 }" m) M出错 Untitled2 (第 20 行)8 ~$ q& j# o8 e7 K- P
[t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
/ z9 \; ^# d" j1 ~* X& V
& E* K6 Q+ Y8 n" F. k }7 c8 O) X* T8 J) o* \
感谢大神们！% C, n k7 x0 [
|
|
/1 
发表于 2021-1-5 10:59
收藏
淘帖
支持！
反对！