w1470852 2024-06-24 18:14 采纳率: 66.7%
浏览 9
已结题

matlab无法将表达式转换为双数组怎么解决?

核心是解一个三自由度微分方程组,解出θ关于时间t的变化图像.设定一个变转速的工况,图中g(t)是功率随时间t变化函数,wm是根据功率求得的角速度.图一是主程序,图二是参数和方程,图三是求解函数,系统提示我无法将表达式转换为双数组,想问一下我这个错在哪里有没有修改办法或者建议,万分感谢

clear;clc;close all;
%初始参数,该参数变化,myfun程序要对应调整
syms t;
g(t)=(1460*9.8+1.04*1460*2.67+0.32*2.25*(10+2.67*t)^2/21.15)*(10+2.67*t)/3600;
wm=(g(t)*(10+2.67*t)*9550/(3600*0.4064/(0.93*14)))*pi*2/60;

zi=4;%内转子极对数
zo=14;%外转子极对数
z1=32;%小齿轮齿数
z2=99;%大齿轮齿数
ms=4.62e-3;
% r2=ms*z1/2;%小齿轮分度圆
% r3=ms*z2/2;%大齿轮分度圆
I=z2/z1;%机械齿轮速比
% G=zo/zi;
% ms1=4.12e-3;
G=(zo+zi)/zi;%磁齿轮速比
%一、计算过程
tmax=10;
% tmax=900*2*pi/wm;% 计算时长,按输入轴的周期数计算
bc=2e-5;
y0=[0;0;0;wm;wm/G;wm/G/I];
% aa=0:1e-4:tmax;
% [t,y]=ode45(@myfun3, aa, y0);
[t,y]=odeRK4sys('myfun3',tmax,bc,y0);

figure;plot(t,y(:,1)/G-y(:,2));
title('时域响应图');
xlabel('时间 (s)');
ylabel('扭转振幅(rad)');%磁齿轮的相对弹性角度
figure;plot(t,y(:,2)/I-y(:,3));
title('时域响应图');
xlabel('时间 (s)');
ylabel('扭转振幅(rad)');%机械齿轮的相对弹性转角;
figure;plot(t,y(:,4),t,y(:,5),t,y(:,6));
title('时域响应图');
xlabel('时间 (s)');
ylabel('转速(rad/s)');
function dydt=myfun3(t,y)
syms t;
g(t)=(1460*9.8+1.04*1460*2.67+0.32*2.25*(10+2.67*t)^2/21.15)*(10+2.67*t)/3600;
wm=(g(t)*(10+2.67*t)*9550/(3600*0.4064/(0.93*14)))*pi*2/60;
zi=4;%内转子极对数
zo=14;%外转子极对数
z1=32;%小齿轮齿数
z2=99;%大齿轮齿数
arf0=20/180*pi;%压力角
ms=4.62e-3;
r2=ms*z1/2;%小齿轮分度圆
r3=ms*z2/2;%大齿轮分度圆
rb2=r2*cos(arf0);%小齿轮基圆
rb3=r3*cos(arf0);%大齿轮基圆
I=z2/z1;%机械齿轮速比
G=(zo+zi)/zi;%磁齿轮速比
Tm=(1460*9.8+1.04*1460*2.67+0.32*2.25*(10+2.67*t)^2/21.15)*0.4064/(0.93*14);%电机输入扭矩
TL=Tm*I*G;%变速器负载
%3 转动惯量
IM=0.081197;%电机转子转动惯量kgm2
I0=0.181197;%磁齿轮低速级转动惯量kgm2
I1=0.013663;%磁齿轮高速级转动惯量kgm2
I2=0.0230596;%机械小齿轮转动惯量
I3=1.9825692;%机械大齿轮转动惯量
IL=116;%整车等效转动惯量

%4 刚度与阻尼
%齿轮的阻尼和刚度
w1=wm*G;%机械小齿轮输入转速
kp=10.51e8;
kb=0.25e8;%两个齿轮的刚度曲线为正弦曲线,
km=kp+kb*sin(w1*t*z1);
sg=0.16;%齿轮阻尼比
cm=2*sg*sqrt(kp*I2*I3/(I2*rb3^2+I3*rb2^2));
 b=1e-4;
 xx=rb2*y(2)-rb3*y(3);
 dxx=rb2*y(5)-rb3*y(6);

  if xx-b>0
      fx=xx-b;
  elseif xx+b<0
      fx=xx+b;
  else
      fx=0;
  end
F23=km*fx+cm*dxx;
T0=1.7*Tm*G;
Tc1=0.1*T0*sin(wm*t);
Tc=0.1*(zi*y(4)-(zi+zo)*y(5));
T21=T0*sin(zi*y(1)-(zi+zo)*y(2))+Tc+Tc1;

y(7)=(Tm-T21/G)/(I0+IM);
y(8)=(T21-r2*F23)/(I2+I1);
y(9)=(-TL+r3*F23)/(I3+IL);

dydt=[y(4)
      y(5)
      y(6)
      y(7)
      y(8)
      y(9)];
end
function [t,y] = odeRK4sys(diffeq,tn,h,y0)
t = (0:h:tn)';         %  Column vector of elements with spacing h
nt = length(t);        %  number of steps
neq = length(y0);      %  number of equations that are simultaneously advanced
y = zeros(nt,neq);     %  Preallocate y for speed
y(1,:) = y0(:)';       %  Assign IC. y0(:) is column, y0(:)' is row vector

h2 = h/2;  h3 = h/3;  h6 = h/6;   %  Avoid repeated evaluation of constants    
k1 = zeros(neq,1);  k2 = k1;      %  Preallocate memory for the Runge-Kutta
k3 = k1;  k4 = k1;  ytemp = k1;   %  coefficients and a temporary vector

%  Outer loop for all steps:  j = time step index;  n = equation number index 
for j=2:nt
   told = t(j-1);   yold = y(j-1,:)';   %  Temp variables, yold is column vector
   k1 = feval(diffeq,told,yold);       %  Slopes at the starting point
   for n=1:neq
     ytemp(n) = yold(n) + h2*k1(n);     %  Estimate all y's at midpoint
   end
   k2 = feval(diffeq,told+h2,ytemp);    %  1st estimate of slopes at midpoint
   for n=1:neq
     ytemp(n) = yold(n) + h2*k2(n);     %  2nd estimate of all y's at midpoint
   end
   k3 = feval(diffeq,told+h2,ytemp);    %  2nd estimate of slopes at midpoint
   for n=1:neq
     ytemp(n) = yold(n) + h*k3(n);      %  Estimate y at end point
   end
   k4 = feval(diffeq,told+h,ytemp);     %  Estimate of slopes at endpoint
   for n=1:neq                          %  Simultaneously advance all y's
     y(j,n) = yold(n) + h6*(k1(n)+k4(n)) + h3*(k2(n)+k3(n));
   end
end

  • 写回答

18条回答 默认 最新

  • 专家-赤兔[在线] 全栈领域优质创作者 2024-06-24 18:14
    关注

    引自免费微信小程序:皆我百晓生

    在您的MATLAB代码中,您尝试使用odeRK4sys函数进行四阶Runge-Kutta(RK4)方法求解。但是,在某些步骤中,您似乎试图将表达式转换为双数组,这可能造成问题。

    为了更好地理解您的问题并为您提供帮助,请允许我解释为什么在某些步骤中使用odeRK4sys可能会导致问题。odeRK4sys是一个通用的求解器,用于执行线性微分方程组。它通过迭代计算来更新状态变量,并且通常会自动处理任何可能存在的边界条件。

    然而,在您的代码中,您尝试将特定的表达式转换为双数组,这可能导致问题。例如,当您尝试将y(4)赋值给变量k1时,这是因为在转换过程中可能丢失了原始变量的引用,从而导致意外的行为。

    为了避免这个问题,您可以考虑使用更传统的数值积分方法,如Euler法或其他步进法,而不是直接调用odeRK4sys。这些方法不需要处理复杂的转换过程,而是专注于逐个步骤地求解微分方程组。

    此外,如果您确实需要使用odeRK4sys来进行更复杂的问题求解,请确保正确地使用它,并确保每次迭代都包括所有必要的步骤。这样可以防止意外的行为发生,同时保持代码的清晰性和可读性。

    另外,如果您能提供更多的上下文信息,比如错误的具体信息、使用的数据类型以及您的具体需求,我可以进一步帮助您找到可能的解决方案。

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论 编辑记录
查看更多回答(17条)

报告相同问题?

问题事件

  • 系统已结题 7月3日
  • 已采纳回答 6月25日
  • 创建了问题 6月24日

悬赏问题

  • ¥15 咨询一下有关于王者荣耀赢藏战绩
  • ¥100 求购一套带接口实现实习自动签到打卡
  • ¥50 MacOS 使用虚拟机安装k8s
  • ¥500 亚马逊 COOKIE我如何才能实现 登录一个亚马逊账户 下发新 COOKIE ..我使用下发新COOKIE 导入ADS 指纹浏览器登录,我把账户密码 修改过后,原来下发新COOKIE 不会失效的方式
  • ¥20 玩游戏gpu和cpu利用率特别低,玩游戏卡顿
  • ¥25 oracle中的正则匹配
  • ¥15 关于#vscode#的问题:把软件卸载不会再出现蓝屏
  • ¥15 vimplus出现的错误
  • ¥15 usb无线网卡转typec口
  • ¥30 怎么使用AVL fire ESE软件自带的优化模式来优化设计Soot和NOx?