找回密码
 注册
关于网站域名变更的通知
查看: 536|回复: 3
打印 上一主题 下一主题

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
我的代码解常微分一直说要返回列向量,有大神帮忙看看嘛
$ I6 F. Y! }' e# O: p" W+ y8 {' ^- h* i' V4 u0 a: J2 p
这是代码:1 F1 s1 \# _2 B. v& ]
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
5 K8 X+ o5 n* L( s( Vswitch flow$ t2 q. a* A& k
case 'creep & recovery'
+ [( s6 {* `- C: ]& }4 x        stress0=arg(1);  q! z% a; B; j6 ~. v4 I% m
        risetime=arg(2);
* d/ {" G/ a- O; A        creeptime=arg(3);
* D- ~" R9 k5 p  W        if t<risetime- r0 U/ v! y1 z& W3 d9 r6 S# d  y
            stressrate=stress0/risetime;+ @& V: r6 S: ?4 z( J3 T" E
            stress=stressrate*t;$ [0 f( s) j3 z. k3 n
        elseif t<creeptime- O9 E% j) i+ x$ l- z
            stressrate=0;! R9 P3 r3 k4 @6 K, w
            stress=stress0;
7 V7 S; I3 v, H+ ?% L        elseif t<creeptime+risetime
4 ^( B% o) [/ Q8 Y5 c$ \$ h& b            stressrate=stress0/risetime;3 b, t5 ~; R* F% n8 N' F
            stress=stress0-stressrate*(t-creeptime);
2 V% K( e2 m' {- c1 a. E* {, X        else
- t. s; @' }5 q4 }1 `, R( l3 \            stressrate=0;3 j" M$ N# i% O. X2 o/ P
            stress=0;
' d! S  z" {4 r        end
; c0 Q" j7 h' S0 F) X8 D        gam1=y(1);! B6 x, l8 T7 f0 N
        gam2=y(2);+ W6 ?0 `% o% x5 F8 B. O
        s(1)=gam2;  s) x$ x( n8 H  W+ U1 S
        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);
. F0 n! v1 s2 N1 e! W- v# O3 _/ P1 A, E
9 x# v5 e. u! A4 J) e' g, p
7 [3 K4 Q0 j" y% v3 `3 C% X5 S0 ~* e. @. [
clear all
& v3 t6 n2 P! C1 Mclc
1 a2 Q$ O; ?; Y. c6 G; o% N0 u# Awarning off$ K' J  `  \- A, h

5 |% P* |- f/ C: veta1=10;  % [Pa.s]2 w' w3 y5 ?% E" |7 s
eta2=20;  % Pa.s]
& f% }  ^7 N/ D8 mG1=2;   % [Pa]
  \) y5 U5 t4 G: C9 j; l$ tG2=5;   %[Pa]
/ H9 P9 C) E- |. t6 i( c& F" ^5 \yieldstress=5;  % [Pa]+ u( ~+ ~. M- K+ |
contactstrain=5;# H' u9 I) B- E2 u' E
; K; |5 o2 ?4 d3 B: f

5 F9 f1 t) v9 qflow='creep & recovery';
0 j+ @0 l, W- vstress0=1:2:10;
; D' F& y5 h2 p7 m5 {% V3 ~) z' v6 Erisetime=0.01;0 w  T5 M2 b! ~# q1 j& E- y5 [: V
creeptime=100;
+ Y9 q! ^+ z3 _+ Z* l9 d/ ]for i=1:length(stress0)
( z2 [9 _* B/ q0 A. C4 L6 V5 E- p% o    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);
! P3 P5 o6 v" c    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on$ t# q, N; \) f  E, |
    a{i}=strcatnum2str(stress0(i),'(Pa)');
( C: D% |7 D- a* n' s1 }5 Kend/ Y7 b) ]1 l& S2 _% f5 K/ X
xlin([1e-4 1000]);   * p% n* K% O8 Y1 F3 T8 F# k* H& n
xlabel('Time(s)');      ylabel('Shear strain(-)');1 c8 z3 z  T' p6 N: \* i
legend(a)/ v4 M+ f4 c/ E( A' p  M6 H

$ v& O# T' T, w, Y/ |1 ^) S5 \3 _$ g( [( B2 G3 ~

& r7 s0 Z& ?; H错误使用 odearguments (第 93 行)' e* z1 j- p& a# Y& s
MYWORK 必须返回列向量。
7 H! _1 T# J: l! t& a6 h9 ?+ m" |6 q  F* {+ \
出错 ode15s (第 152 行)
7 q3 {7 w/ X6 X    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);; R3 W/ n* P8 n  s! z
( F, v2 Q! j  d; f
出错 Untitled2 (第 20 行)
" I# s6 I6 l3 T- o! n* j# D3 T/ ]    [t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);; v/ r1 T8 p* A, ~4 A/ M
) [% T' i- Y  a% b
; {( t! a! x5 U- y
感谢大神们!
+ `$ A, g# _- ^! D  r# b

该用户从未签到

2#
发表于 2021-1-5 13:09 | 只看该作者
clear all
1 T0 ]4 H, X5 K' d# T% m" lclc
, P' w& {8 M" u5 V& _warning off! Y4 J3 X: r& _. H

& ~1 L% H& C$ U0 O/ n4 A0 [, zeta1=10;  % [Pa.s]
) [9 |( i/ d+ Aeta2=20;  % Pa.s]! N0 f! [% `+ z' [
G1=2;   % [Pa]
) V4 m! H  o! o! M( S6 G2 `, lG2=5;   %[Pa]% a0 _' G+ l& n" `, [3 d! A
yieldstress=5;  % [Pa]
: `1 C" s9 T1 Y. @: L& tcontactstrain=5;* k, }+ u! c- d+ ~) \% s# J
% t5 K; {" c- J( e

& s  N) O2 _7 u& e, Mflow='creep & recovery';
/ ^( t  V( y' K" Xstress0=1:2:10;
' r0 Z# {7 |4 X5 Irisetime=0.01;
% E; T+ ?7 c, a$ @" W. n# N# qcreeptime=100;
5 B5 x( w- V; x5 dfor i=1:length(stress0)
5 i. o( l4 Q% D& ?. I. c8 W    [t,strain]=ode15s(@mywork,[0 10],[0 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);/ a6 `+ S: n8 M& H' r1 c
    subplot(2,3,3);     semilogx(t,strain(:,1));        hold on
: m3 ^1 j+ Q! R! ]2 B    a{i}=strcat(num2str(stress0(i)),'(Pa)');3 t. Z2 u3 O' M# F
end
" l0 Q5 |1 [$ E, u2 C/ X, Nxlim([1e-4 1000]);" A# \& v# D1 m  B. p
xlabel('Time(s)');      ylabel('Shear strain(-)');+ ]1 ]! I7 U; d2 I, Z
legend(a)$ V0 ~! S3 k9 g; S8 Q# B
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
  v& ]9 u6 `' t7 q( m$ b$ cswitch flow
$ ~& w5 Y) p1 V$ H# Y0 H    case 'creep & recovery'
9 K! S! D% j. Y        stress0=arg(1);4 t  l! z: J+ C) X9 i& E
        risetime=arg(2);
4 [# s1 {6 X' O  ]6 z; I9 m        creeptime=arg(3);
7 v4 T! V6 ~6 J/ @2 e0 p        if t<risetime
9 X; O  {4 r) B( }            stressrate=stress0/risetime;: Q  |/ A4 ]$ V: E
            stress=stressrate*t;% D- |* l) H, F
        elseif t<creeptime
) ]6 m9 |- O2 B* M; D8 v            stressrate=0;
% r( A7 D2 B3 l. W/ n# @            stress=stress0;
! y- X  n- ?# k        elseif t<creeptime+risetime% |6 R" q' s' T4 X6 L4 h( N
            stressrate=stress0/risetime;
( n3 q8 B+ ?& w; F8 t            stress=stress0-stressrate*(t-creeptime);
( k- p3 ?0 A  l" F/ v        else
- D0 z3 ^) l: R5 |: |( S            stressrate=0;
/ R, S4 V0 c  X, m5 I) `0 Z            stress=0;9 L% I- b$ y+ y1 P( ?2 n9 {8 w4 {+ K
        end; \0 w( K! S8 [2 T; K+ i: L
        gam1=y(1);
1 U$ c, x+ E2 m0 @4 \1 C( T        gam2=y(2);6 R' ~8 i7 k" C* x* H7 K* R" J! Z
        s(1)=gam2;
- t/ H$ U) y! @* Z  d+ w% F& x        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 j- @3 y* j8 E* g, R$ Z% j
        s=s';: X: y) r" g& U+ I6 I5 J) T
end
: R$ t0 D4 h$ ^6 E! T, [/ Mend

该用户从未签到

4#
发表于 2021-1-5 15:43 | 只看该作者
来学习一下
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-10-31 19:46 , Processed in 0.125000 second(s), 23 queries , Gzip On.

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

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

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