[已答复] 大家好,请问MATLAB程序运行出不来图,是什么原因,救急谢谢啦!

[复制链接]
KellyGX001 发表于 2022-5-13 17:50:43
附件中有源程序


function dxx=diyicihanshu(~,x)
global  F  W

x1=x(1);%主振子m1位移
x2=x(2); %速度
x3=x(3);%吸振器m2位移
x4=x(4);
x5=x(5);%m3的位移
x6=x(6);
x7=x(7);%m4的位移
x8=x(8);
x9=x(9);


k1=1500;
k2=300;
k3=300;
km=5000;
cc1=0.4;
cc2=0.1;
m1=3;
m2=0.3;
m3=0.05;
m4=0.05;
br=1.26;%剩余磁通密度
A2=0.01.^2*pi;%活动圆形面面积
A1=0.01.^2*pi;
% a2=0.05;%磁铁之间的水平间隙
% a1=0.05;
r2=0.01;%圆柱形磁铁的半径
r1=0.01;
p=0.05;%磁铁厚度
u=1.256*1e-6;%磁导率

dx1=x2;
dx2=-(k1/m1)*x1-(cc1/m1)*x2-(cc2/m1)*(x2-x4)-km*(x1-x3).^3-k2*(x1-x5)+k3*(x1-x7)+(F/m1)*cos(W*x9);
% dx2=-(k1/m1)*x1-(cc1/m1)*x2-(cc2/m1)*(x2-x4)-km*(x1-x3).^3-k2*(x1-x5)+k3*(x1-x7)+(F/m1)*cos(W*x9)-(br.^2*A2.^2*(t+r2).^2./(m1*pi*u*t.^2))*(1./(a2+x3).^2+1./(a2+x3+2*t).^2-2./(a2+x3+t).^2)+(br.^2*A1.^2*(t+r1).^2./(m1*pi*u*t.^2))*(1./(a1-x3).^2+1./(a1-x3+2*t).^2-2./(a1-x3+t).^2);
dx3=x4;
% dx4=(cc2/m2)*(x2-x4)+km*(x1-x3).^3;
dx4=(cc2/m2)*(x2-x4)+km*(x1-x3).^3+(br.^2*A2.^2*(p+r2).^2./(m1*pi*u*p.^2))*(1./(x3-x5).^2+1./(x3-x5+2*p).^2-2./(x3-x5+p).^2)-(br.^2*A1.^2*(p+r1).^2./(m1*pi*u*p.^2))*(1./(x7-x3).^2+1./(x7-x3+2*p).^2-2./(x7-x3+p).^2);
dx5=x6;
% dx6=(k2/m3)*(x1-x5);
dx6=(k2/m3)*(x1-x5)-(br.^2*A2.^2*(p+r2).^2./(m1*pi*u*p.^2))*(1./(x3-x5).^2+1./(x3-x5+2*p).^2-2./(x3-x5+p).^2);
dx7=x8;
% dx8=(k3/m4)*(x1-x7);
dx8=(k3/m4)*(x1-x7)+(br.^2*A1.^2*(p+r1).^2./(m1*pi*u*p.^2))*(1./(x7-x3).^2+1./(x7-x3+2*p).^2-2./(x7-x3+p).^2);
dx9=1;
dxx=[dx1,dx2,dx3,dx4,dx5,dx6,dx7,dx8,dx9]';





clc;
clear all;

opts = odeset('RelTol',1e-5,'AbsTol',1e-3,'MaxStep',0.1);

global    F W

tspan=0.1:0.01:60;
y0=[0 0 0 0 0 0 0 0 0]';
F=10;
% y0=[0.522465733302094        2.18583198169893        1.60287314982309        0.759192863254238        0];
% y0=[0.01 0.01 0.01 0.1 0];
% W=13;

%龙格库塔法相当于在tspan区间内,用一个tspan(i),计算出一个 y(i,:),用其他方法,设定I为2s一变,在0.1(弱非)与0.2(强非)之间;
W=10;
% tic;
[t,y] = ode45('diyicihanshu', tspan, y0,opts);

% time=toc;
% hold on
% plot((5-t/40000)*6.28,y(:,1),'r','LineWidth',1);
plot(t,y(:,1),'k','LineWidth',1);
% plot(t,y(:,1)-y(:,3),'k','LineWidth',1);
% plot(t,y(:,1)-y(:,3),'k','LineWidth',1);
set(gca,'FontSize',18,'FontWeight','bold')
xlabel('\itt/\rms','FontSize',30,'FontWeight','bold','FontName','Times New Roman');
ylabel('\itw/\rmm','FontSize',30,'FontWeight','bold','FontName','Times New Roman')

微信截图_20220513162939.png

diyicihanshu.m

1.44 KB, 下载次数: 1

自定义函数

diyici.m

879 Bytes, 下载次数: 1

1 条回复


wjb986555360 发表于 5 天前
求的y值都是NAN值,你的方程可能写的有问题
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表