幅频响应曲线的ＭＡＴＬＡＢ画法

mechanic05

各位同人暑期过的可好？
我想请教关于幅频响应曲线的ＭＡＴＬＡＢ作图方法．我的幅频响应方程为附件．我是直接用ＭＡＰＬＥ或者ＭＡＴＬＡＢ中的ＰＬＯＴ语句画的图，很不光滑．看外文文献上人家说根据幅频响应方程用ＭＡＴＬＡＢ作图就得到有跳跃和多值的幅频响应曲线．　　
麻烦前辈能给指点！

mechanic05

:'(没人理我．我在尝试新的方法，但愿成功，到时一定和大家分享，讨论．:victory::handshake

lzjtu

解出a和a1的关系式，直接用plot就可以画出来了

mechanic05

谢谢！可是画出来的很不光滑．也不知道韦什麽？请指教。:@)

lzjtu

把步长取小点试试

21172485

一般来说这是高次代数方程,振幅无法显式表达,但是你的siagma1(你的方程里表示频率的)是2次方,是可以求解的,那么求解以后
用MATLAB求解.总之就是哪个次数低求那个.另外用小步长其实也不是一个好的办法,因为那样也无法太光滑,并且耗时太多
直接求解比较简单,这也是我硕士论文期间想出的办法
clc
clear
u1=-0.5;
u2=-0.1;
t1=0.2;
t2=30;
q1=50000;
a=0:0.01:4.075
p=u1.*a+u2*a.^2-t1*a.^5
c1=sqrt(q1-p.^2)+t2*a.^3
c2=-sqrt(q1-p.^2)+t2*a.^3
plot(c1,a)
hold on
plot(c2,a)

mechanic05

:@) 非常感谢前辈！！
这两天忙别的事情都没顾上上论坛．我也画了图，趋势有些像，可太欠光滑．我得仔细研究下前辈的程序．向前辈致敬！！

mechanic05

大家好！我把这个幅频响应图用ＭＡＰＬＥ重新又画了一下，把步长取的小了一些，效果还可以．ＭＡＴＬＡＢ和ＭＡＰＬＥ是一样的，取小步长效果应该不错．现在把程序和图贴出来供大家参考．希望共同提高．
fupin222:=(-A*omega0*(omega-omega0)+3*omega0^2*epsilon*A^3/8)^2+zeta^2*A^2*omega0^4=F^2/4;
                                                      2          3 2
  fupin222 := (-A omega0 (omega - omega0) + 3/8 omega0  epsilon A )
               2  2       4        2
         + zeta  A  omega0  = 1/4 F
> epsilon:=0.1: F:=0.1: zeta:=0.01: omega0:=1:
> with(plots,implicitplot):
> implicitplot(fupin222,omega=0..4,A=0..6,numpoints=10000);

无水1324

不错的解决方法

无水1324

这个图还是不连贯哈

21172485

你使用了隐函数绘图函数,结果是对的
但是效果不是很好.呵呵.缩小步长还是许多点的集合而已

mechanic05

前辈好！
就我的问题我首先用隐函数画图，如图一；接着又用前辈的方法画图，如图二．看起来前辈的方法还是颇有一些优点的，比如可以对不同的部分用不同的图线．

mechanic05

mechanic05

我把上面这附图的程序附上。麻烦前辈给修改一下，就是把稳定区与不稳定区域能用不同的线型画出．我只能作到上面那样．跳跃部分是实际不存在的运动．应该和前后两部分的稳定区域区别开．试验半天也不行．麻烦前辈给指点哈！谢谢！
clc
clear
u1=0.3;
u2=0.25;
w0=1;
gv=0.3;
a1=4;a2=0.5;a3=0.2;
b1=1;
t1=5*gv./8;
t2=(3*a1-3*a2+a3)./8;
q1=b1./2;
a=0.1:0.001:1.53;
p=u1+u2*a-t1*a.^4;
m=q1./a;
k=sqrt(m.^2-p.^2);
f1=(t2*a.^2-k);
f2=(t2*a.^2+k);
plot(f1,a,'.');
hold on
plot(f2,a,':');

mechanic05

不好意思，ＭＡＴＬＡＢ不熟悉．
为什摸这样改不行呢？
clear
u1=0.3;
u2=0.25;
w0=1;
gv=0.3;
a1=4;a2=0.5;a3=0.2;
b1=1;
t1=5*gv./8;
t2=(3*a1-3*a2+a3)./8;
q1=b1./2;
a=0.1:0.001:1.53;
p=u1+u2*a-t1*a.^4;
m=q1./a;
k=sqrt(m.^2-p.^2);
f1=(t2*a.^2-k);
f2=(t2*a.^2+k);
plot(f1,a,'.');
hold on
a=[0.1:0.001:0.17];
plot(f1,a,':');
hold on
a=[0.17:0.001:1.53];
plot(f1,a,'.')
请指点！