本站微信 收藏本站 设为首页
 找回密码
 注册
关于网站域名变更的通知
EDA365电子论坛网»电子社区 › 研发技术版块 › MATLAB论坛 › matlab解常微分一直报错返回列向量
返回列表 发新帖
查看: 547|回复: 3
打印 上一主题 下一主题

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

[复制链接]
cichishia

该用户从未签到
跳转到指定楼层
1#
发表于 2021-1-5 10:59 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录！

您需要 登录 才可以下载或查看，没有帐号？注册

x
我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛
7 l; p; q$ C3 V! ]$ c9 A( k* @/ d7 g" _6 J3 U0 V
这是代码：( X2 b2 b6 u* l* i. t
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
( {' `$ k4 M' S% |- s: Iswitch flow4 \/ G* `6 o6 N! X5 L' A
case 'creep & recovery'
# i1 s. l+ k* ?5 d        stress0=arg(1);
0 O' L1 c1 P; m# [9 S7 S% I        risetime=arg(2);/ I9 z: T9 `  G
        creeptime=arg(3);
6 z0 Z9 ^0 m$ p+ c        if t<risetime( O; ]0 L3 h. p) N4 _  v, O
            stressrate=stress0/risetime;; Z* W$ n2 s$ J/ k2 T
            stress=stressrate*t;4 P* c, l9 X1 R( z/ u) E( S
        elseif t<creeptime& R7 Q; k' i) C* G1 a& Y
            stressrate=0;
4 O% J. K2 F6 O! p6 ^' a            stress=stress0;  M) f$ ?/ H3 _. k6 c. M: P
        elseif t<creeptime+risetime
1 @) e! d, M2 o) B! M% y$ u: z  o            stressrate=stress0/risetime;
7 t& D3 `! y7 Z            stress=stress0-stressrate*(t-creeptime);( t. `0 F4 m8 `
        else- s: R- B  u# K7 e" s
            stressrate=0;! C1 X4 b8 c5 ~* W3 Q) u" }( D$ y
            stress=0;# M9 V+ T* z' O+ }2 q1 Q
        end" G9 H& N5 i4 x/ V9 e
        gam1=y(1);8 D0 @/ m% |; W: R+ v
        gam2=y(2);2 Q5 q+ n# @4 L: V5 r' I
        s(1)=gam2;
  T9 K" n0 E) K        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);9 C$ U$ J- @3 _4 k6 r2 D
1 T- I' ]  @8 y  ~) K
$ x) X9 |0 u, r! e; Q, h6 U8 \
  n0 F# l# G! k: k  H: y& D
clear all
; p  l/ E6 J+ u4 U% u1 ^" l. gclc7 h: B; Z6 u- s; @  G
warning off
- m( W9 L) }* J9 Q# p3 N7 N0 f8 U: d; i6 s# \
eta1=10;  % [Pa.s]" n" b- X% _! K1 j8 C- v
eta2=20;  % Pa.s]
7 o# }# b% z, z. I# Q) e: h: gG1=2;   % [Pa]
% \# l1 E2 H& z; d2 {G2=5;   %[Pa]
' e$ l3 g1 N9 R1 v+ yyieldstress=5;  % [Pa], z8 M/ x* F/ _0 a, _3 k
contactstrain=5;
& O& C7 P* y4 n7 A0 F4 O, J2 J) m( @
6 {: {/ p" W0 ?% W6 E! _& P5 m. H! [$ n* v1 [& t( A- Z1 m  P& a
flow='creep & recovery';
- m3 c/ X7 ?. }1 o" Z6 mstress0=1:2:10;+ q0 ]" }% y5 a( Q8 H# Y
risetime=0.01;
: {- r7 t5 c1 x* U( m- wcreeptime=100;! F$ [) k/ c( m6 l- G0 [5 ]2 E  Z
for i=1:length(stress0)
' ^, D5 B5 q7 d    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);0 ?* |5 T% J) _+ ~  }
    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
& l- |: V$ Y% X% k5 i, S    a{i}=strcatnum2str(stress0(i),'(Pa)');* c- Q! h& y# p9 g' G$ H! S
end
& w& x7 `* K$ L7 U1 [2 Gxlin([1e-4 1000]);   , b* b8 A) F2 @+ O
xlabel('Time(s)');      ylabel('Shear strain(-)');
5 y+ Z, ~! p" L2 i; {legend(a)3 }4 q2 |6 \+ O4 H  W

8 ~% H8 Q+ V) \; ]' c2 T! D7 r  X7 T- x9 |
0 E7 R7 W( `: p8 |1 `$ Z& j+ a
错误使用 odearguments (第 93 行), f2 D+ r. r. f- l( t
MYWORK 必须返回列向量。
+ ?; \5 @5 J5 t6 [( m* d( s: p0 a4 y- x
出错 ode15s (第 152 行)
9 R' f' p# G" o1 ^2 v    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);7 [$ y' b8 M- |! G6 Q( I! N

5 }" @  Y2 c) }5 R出错 Untitled2 (第 20 行)
+ u' o3 l6 }9 _) W; s    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);6 z* T! m5 E; ]8 |% }! i; E
7 @- }/ }, i1 }8 e9 a- X+ T* L
0 m( K; A  Y1 d3 n
感谢大神们！
8 l, W+ R1 F& y9 t2 S" t0 n' m

下载资料威望不够？点击查看获取威望的N种方法>>
回复

举报

EDA365微信号及QQ群号！
llbnmo

该用户从未签到
2#
发表于 2021-1-5 13:09 | 只看该作者
clear all
1 ]0 {3 P7 P8 D. ^clc
! |4 E8 P: p$ O6 ewarning off% v* `8 c% N& Y1 E3 f  v

# ^( {7 @0 @5 Z7 F" m, C# O7 ~eta1=10;  % [Pa.s]" U5 H" i) ]& ?' p8 f1 o
eta2=20;  % Pa.s]( D  v3 h' L( R, x* R$ b- `
G1=2;   % [Pa]
7 Q* M# I/ T* u9 o' u* S9 \G2=5;   %[Pa]! ?) G. Z- n+ u- f# `5 j
yieldstress=5;  % [Pa], n& C7 ^' H- X; n' F1 e6 A
contactstrain=5;
6 p- h9 s0 W) ]8 h( u
; e) c' G. m$ ~# e7 o; q4 R* k9 O4 g7 Y+ ]2 ]
flow='creep & recovery';
; K% R+ N# G' _6 i2 {' {' dstress0=1:2:10;- T* D. C  R8 u. }2 t/ Q
risetime=0.01;
. L4 w$ e" n7 b2 P% H! [creeptime=100;) |* B$ V& t3 z, [8 |# D
for i=1:length(stress0)9 Z) i; e% P3 t' I- e1 e4 u( E
    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);) b' d5 R' H4 X4 X0 S( }
    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on$ P+ Z. o0 s6 D* n) {# V1 _
    a{i}=strcat(num2str(stress0(i)),'(Pa)');
+ ~3 k0 Z# l1 R: K5 z+ I; pend, C: q6 n6 d2 v$ m" _1 z
xlim([1e-4 1000]);
3 ?2 o& n; g/ Y% Zxlabel('Time(s)');      ylabel('Shear strain(-)');
7 Q# ]& r" W7 Z& T6 _& C6 q- j; glegend(a)
+ ?! b  N! s" w) o/ `function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
2 A1 f% I2 U' k7 t$ tswitch flow3 ?% H7 D2 h" w& G
    case 'creep & recovery'& r; E3 b5 I( J# P" H" F8 S1 A
        stress0=arg(1);4 {( s6 h* J0 X# `
        risetime=arg(2);6 i# k# [$ p2 G) I- E) c4 w/ Z7 L
        creeptime=arg(3);& A' g2 U3 H0 K5 l
        if t<risetime* F; V. u# q, i% {& z( a
            stressrate=stress0/risetime;# L! c8 B1 l: p
            stress=stressrate*t;
" C2 \/ I. r9 u6 e        elseif t<creeptime9 ]5 d1 S! U8 s& `  i+ [3 _' V
            stressrate=0;$ p1 S# w$ ]( M6 w0 l
            stress=stress0;
% z0 g' g( |7 r; e8 o        elseif t<creeptime+risetime
+ n% X% \1 {3 E% X            stressrate=stress0/risetime;
/ B2 M4 B# m4 b* q! L! S& c) B            stress=stress0-stressrate*(t-creeptime);& {' h8 S7 H8 ~; J" \
        else
) U, K  `$ v8 y6 O  q) R) G1 G            stressrate=0;
& r' q6 w" @! L- V$ E8 z9 x            stress=0;
% q; d- I' v8 g        end" C9 ^4 U( ~7 L* t# i6 S
        gam1=y(1);
/ J! k. ^0 r8 |( S2 H        gam2=y(2);
" |5 T" T' @( d2 m+ }, S        s(1)=gam2;
' h/ W) l) b/ |5 t        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);3 V- _3 x) S6 L9 _% _+ Y2 u
        s=s';  ^; u6 m* F1 W" V7 A
end
6 M/ D9 N1 K) |3 wend

下载资料威望不够？点击查看获取威望的N种方法>>
回复 支持 反对

举报

IBB-EUT

该用户从未签到
4#
发表于 2021-1-5 15:43 | 只看该作者
来学习一下

下载资料威望不够？点击查看获取威望的N种方法>>
回复 支持 反对

举报

返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-11-24 16:18 , Processed in 0.187500 second(s), 24 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表