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

cichishia · 发表于 2021-1-5 10:59

EDA365欢迎您登录！

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

x

我的代码解常微分一直说要返回列向量，有大神帮忙看看嘛

这是代码：
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
case 'creep & recovery'
      stress0=arg(1);
      risetime=arg(2);
      creeptime=arg(3);
      if t<risetime
         stressrate=stress0/risetime;
         stress=stressrate*t;
      elseif t<creeptime
         stressrate=0;
         stress=stress0;
      elseif t<creeptime+risetime
         stressrate=stress0/risetime;
         stress=stress0-stressrate*(t-creeptime);
      else
         stressrate=0;
         stress=0;
      end
      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);

clear all
clc
warning off

eta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
G1=2; % [Pa]
G2=5; %[Pa]
yieldstress=5;  % [Pa]
contactstrain=5;

flow='creep & recovery';
stress0=1:2:10;
risetime=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]);
subplot(2,3,3);    semilogx(t,strain(:,1));       hold on
a{i}=strcatnum2str(stress0(i),'(Pa)');
end
xlin([1e-4 1000]);
xlabel('Time(s)');    ylabel('Shear strain(-)');
legend(a)

错误使用 odearguments (第 93 行)
MYWORK 必须返回列向量。

出错 ode15s (第 152 行)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

出错 Untitled2 (第 20 行)
[t,strain]=ode15s(@mywork,[0 10],[1 0],[],eta1,eta2,G1,yieldstress,contactstrain,flow,[stress0(i) risetime creeptime]);

感谢大神们！

llbnmo · 发表于 2021-1-5 13:09

clear all
clc
warning off

eta1=10;  % [Pa.s]
eta2=20;  % Pa.s]
G1=2; % [Pa]
G2=5; %[Pa]
yieldstress=5;  % [Pa]
contactstrain=5;

flow='creep & recovery';
stress0=1:2:10;
risetime=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]);
subplot(2,3,3);    semilogx(t,strain(:,1));       hold on
a{i}=strcat(num2str(stress0(i)),'(Pa)');
end
xlim([1e-4 1000]);
xlabel('Time(s)');    ylabel('Shear strain(-)');
legend(a)
function s=mywork(t,y,eta1,eta2,G1,yieldstress,contactstrain,flow,arg)
switch flow
case 'creep & recovery'
      stress0=arg(1);
      risetime=arg(2);
      creeptime=arg(3);
      if t<risetime
         stressrate=stress0/risetime;
         stress=stressrate*t;
      elseif t<creeptime
         stressrate=0;
         stress=stress0;
      elseif t<creeptime+risetime
         stressrate=stress0/risetime;
         stress=stress0-stressrate*(t-creeptime);
      else
         stressrate=0;
         stress=0;
      end
      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);
      s=s';
end
end

小白的白 · 发表于 2021-1-5 13:17

楼上正解！

IBB-EUT · 发表于 2021-1-5 15:43

来学习一下

帐号		自动登录	找回密码
密码			注册

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

EDA365欢迎您登录！

浏览过的版块

推荐内容 /1