查看: 322|回复: 3|关注: 0

[已解决] 请问一下matlab绘制动态曲线数据精度错误的问题

[复制链接]

新手

7 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
  • 关注者: 1
本帖最后由 328002469 于 2020-7-31 18:33 编辑

请直接拉到最下面代码查看。所要绘制的就是r1,r2,r3三个解关于时间t的动态曲线,按照help里的命令改成r1后就报错:参数 Y 的类型无效。类型必须为双精度。请问这种情况该如何解决。刚接触matlab没多久,希望大家帮忙解答一下谢谢。


clc;
clear all;
M=[9 0 0;0 1 0;0 0 3];%M=[m1 0 0;0 m2 0;0 0 m3]
C=[3 3 0;3 3.3 -0.3;0 -0.3 0.3];%C=[c1 c1 0;c1 c1+c2 -c2;0 -c2 c2]
K=[27 27 0;27 30 -3;0 -3 3];%K=[k1 k1 0;k1 k1+k2 -k2;0 -k2 k2]
x0=[1;0;0];
dx0=[0;0;0];%条件矩阵
xa=0;
xb=0;
xc=0;%定义数值积分求解变量
Ms=(chol(M))^(-1);%求M的^(-1/2)
Cw=Ms*C*Ms;
Kw=Ms*K*Ms;
[P,t1]=eig(Kw);%求Kw的特征向量并赋予给P,用来对角化Kw
PT=inv(P);%求P的转置矩阵
Kh=P\Kw*P;
Ch=P\Cw*P;%分别将Kw和Cw对角化
X=P\Ms;
syms t;
F1=2*sin(t);
F2=3*cos(t);
F3=3*sin(t);
F=[F1;F2;F3];
f=X*F; %函数矩阵
S=Ms*P;%求坐标变换矩阵
Sd=S^(-1);
r0=Sd*x0;
dr0=Sd*dx0;%求在解耦条件下的初始条件
Ch11=Ch(1,1);
Ch22=Ch(2,2);
Ch33=Ch(3,3);
Kh11=Kh(1,1);
Kh22=Kh(2,2);
Kh33=Kh(3,3);%取出矩阵某行某列的数
f1=f(1,1);
f2=f(2,1);
f3=f(3,1);
f11=subs(f1,t,0);
f22=subs(f2,t,0);
f33=subs(f3,t,0);
syms r1(t) r2(t) r3(t);%定义表示微分方程结果的符号变量
eqn1=diff(r1,2)+Ch11*diff(r1)+Kh11*r1==f1;
Dr1=diff(r1);
cond1=[r1(0)==f11,Dr1(0)==0];
eqn2=diff(r2,2)+Ch22*diff(r2)+Kh22*r2==f2;
Dr2=diff(r2);
cond2=[r2(0)==f22,Dr2(0)==0];
eqn3=diff(r3,2)+Ch33*diff(r3)+Kh33*r3==f3;
Dr3=diff(r3);
cond3=[r3(0)==f33,Dr3(0)==0];
r1=dsolve(eqn1,cond1);%求解微分方程
r2=dsolve(eqn2,cond2);%r1,r2,r3分别为三个质量块的振动响应
r3=dsolve(eqn3,cond3);


r1=simplify(r1);
r2=simplify(r2);
r3=simplify(r3);



figure;
subplot(1,2,1);
C1=0;D1=0;G1=0;
C1=eval(r1);%用eval函数引入r1的值
h = animatedline;
axis([0,30,-15,10])

t = linspace(0,30,1000);
y = C1;
for k = 1:length(t)
    addpoints(h,t(k),y(k));
    drawnow
end


回复主题 已获打赏: 0 积分

举报

论坛优秀回答者

9

主题

1635

帖子

343

最佳答案
  • 关注者: 81
发表于 2020-7-31 19:28:00 | 显示全部楼层 |此回复为最佳答案
仅供参考
  1. clc;
  2. clear all;
  3. M=[9 0 0;0 1 0;0 0 3];%M=[m1 0 0;0 m2 0;0 0 m3]
  4. C=[3 3 0;3 3.3 -0.3;0 -0.3 0.3];%C=[c1 c1 0;c1 c1+c2 -c2;0 -c2 c2]
  5. K=[27 27 0;27 30 -3;0 -3 3];%K=[k1 k1 0;k1 k1+k2 -k2;0 -k2 k2]
  6. x0=[1;0;0];
  7. dx0=[0;0;0];%条件矩阵
  8. xa=0;
  9. xb=0;
  10. xc=0;%定义数值积分求解变量
  11. Ms=(chol(M))^(-1);%求M的^(-1/2)
  12. Cw=Ms*C*Ms;
  13. Kw=Ms*K*Ms;
  14. [P,t1]=eig(Kw);%求Kw的特征向量并赋予给P,用来对角化Kw
  15. PT=inv(P);%求P的转置矩阵
  16. Kh=P\Kw*P;
  17. Ch=P\Cw*P;%分别将Kw和Cw对角化
  18. X=P\Ms;
  19. syms t;
  20. F1=2*sin(t);
  21. F2=3*cos(t);
  22. F3=3*sin(t);
  23. F=[F1;F2;F3];
  24. f=X*F; %函数矩阵
  25. S=Ms*P;%求坐标变换矩阵
  26. Sd=S^(-1);
  27. r0=Sd*x0;
  28. dr0=Sd*dx0;%求在解耦条件下的初始条件
  29. Ch11=Ch(1,1);
  30. Ch22=Ch(2,2);
  31. Ch33=Ch(3,3);
  32. Kh11=Kh(1,1);
  33. Kh22=Kh(2,2);
  34. Kh33=Kh(3,3);%取出矩阵某行某列的数
  35. f1=f(1,1);
  36. f2=f(2,1);
  37. f3=f(3,1);
  38. f11=subs(f1,t,0);
  39. f22=subs(f2,t,0);
  40. f33=subs(f3,t,0);
  41. syms r1(t) r2(t) r3(t);%定义表示微分方程结果的符号变量
  42. eqn1=diff(r1,2)+Ch11*diff(r1)+Kh11*r1==f1;
  43. Dr1=diff(r1);
  44. cond1=[r1(0)==f11,Dr1(0)==0];
  45. eqn2=diff(r2,2)+Ch22*diff(r2)+Kh22*r2==f2;
  46. Dr2=diff(r2);
  47. cond2=[r2(0)==f22,Dr2(0)==0];
  48. eqn3=diff(r3,2)+Ch33*diff(r3)+Kh33*r3==f3;
  49. Dr3=diff(r3);
  50. cond3=[r3(0)==f33,Dr3(0)==0];
  51. r1=dsolve(eqn1,cond1);%求解微分方程
  52. r2=dsolve(eqn2,cond2);%r1,r2,r3分别为三个质量块的振动响应
  53. r3=dsolve(eqn3,cond3);


  54. r1=simplify(r1);
  55. r2=simplify(r2);
  56. r3=simplify(r3);



  57. figure;
  58. subplot(1,2,1);
  59. C1=0;D1=0;G1=0;
  60. C1=eval(r1);%用eval函数引入r1的值
  61. h = animatedline;
  62. axis([0,30,-15,10])

  63. t1 = linspace(0,30,1000);
  64. C1=double(subs(C1,t,t1));
  65. y = C1;
  66. for k = 1:length(t1)
  67.     addpoints(h,t1(k),y(k));
  68.     drawnow
  69. end
复制代码
回复此楼 已获打赏: 0 积分

举报

新手

7 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
  • 关注者: 1
 楼主| 发表于 2020-7-31 19:57:34 | 显示全部楼层

谢谢,问题已解决
回复此楼 已获打赏: 0 积分

举报

新手

7 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
  • 关注者: 1
 楼主| 发表于 2020-7-31 20:13:47 | 显示全部楼层

还想再问一下您,我用subplot(1,2,1)和subplot(1,2,2)创立了两幅动态图,执行绘制命令的时候是左面画完再画右面,如何让左右两幅图同时进行绘制呢。

clc;
clear all;
M=[9 0 0;0 1 0;0 0 3];%M=[m1 0 0;0 m2 0;0 0 m3]
C=[3 3 0;3 3.3 -0.3;0 -0.3 0.3];%C=[c1 c1 0;c1 c1+c2 -c2;0 -c2 c2]
K=[27 27 0;27 30 -3;0 -3 3];%K=[k1 k1 0;k1 k1+k2 -k2;0 -k2 k2]
x0=[1;0;0];
dx0=[0;0;0];%条件矩阵
xa=0;
xb=0;
xc=0;%定义数值积分求解变量
Ms=(chol(M))^(-1);%求M的^(-1/2)
Cw=Ms*C*Ms;
Kw=Ms*K*Ms;
[P,t1]=eig(Kw);%求Kw的特征向量并赋予给P,用来对角化Kw
PT=inv(P);%求P的转置矩阵
Kh=P\Kw*P;
Ch=P\Cw*P;%分别将Kw和Cw对角化
X=P\Ms;
syms t;
F1=2*sin(t);
F2=3*cos(t);
F3=3*sin(t);
F=[F1;F2;F3];
f=X*F; %函数矩阵
S=Ms*P;%求坐标变换矩阵
Sd=S^(-1);
r0=Sd*x0;
dr0=Sd*dx0;%求在解耦条件下的初始条件
Ch11=Ch(1,1);
Ch22=Ch(2,2);
Ch33=Ch(3,3);
Kh11=Kh(1,1);
Kh22=Kh(2,2);
Kh33=Kh(3,3);%取出矩阵某行某列的数
f1=f(1,1);
f2=f(2,1);
f3=f(3,1);
f11=subs(f1,t,0);
f22=subs(f2,t,0);
f33=subs(f3,t,0);
syms r1(t) r2(t) r3(t);%定义表示微分方程结果的符号变量
eqn1=diff(r1,2)+Ch11*diff(r1)+Kh11*r1==f1;
Dr1=diff(r1);
cond1=[r1(0)==f11,Dr1(0)==0];
eqn2=diff(r2,2)+Ch22*diff(r2)+Kh22*r2==f2;
Dr2=diff(r2);
cond2=[r2(0)==f22,Dr2(0)==0];
eqn3=diff(r3,2)+Ch33*diff(r3)+Kh33*r3==f3;
Dr3=diff(r3);
cond3=[r3(0)==f33,Dr3(0)==0];
r1=dsolve(eqn1,cond1);%求解微分方程
r2=dsolve(eqn2,cond2);%r1,r2,r3分别为三个质量块的振动响应
r3=dsolve(eqn3,cond3);


r1=simplify(r1);
r2=simplify(r2);
r3=simplify(r3);



figure;
subplot(1,2,1);
C1=0;D1=0;G1=0;
C1=eval(r1);%用eval函数引入r1的值
D1=eval(r2);
G1=eval(r3);

h1=animatedline('Color','b','LineStyle','-','LineWidth',2);
h2=animatedline('Color','r','LineStyle','-.','LineWidth',2);
h3=animatedline('Color','g','LineStyle','--','LineWidth',2);
axis([0,30,-15,10])
set(gca,'XTick',0:30/5:30); %设定x轴坐标刻度
set(gca,'YTick',-15:25/5:10); %设定y轴坐标刻度
xlabel('时间(T/s)');
rectangle('position',[5,9,1,0.1],'FaceColor','b');
text(0,9.1,'质量块1振动响应','Color','black','fontsize',10);
rectangle('position',[5,8.5,1,0.1],'FaceColor','r');
text(0,8.6,'质量块2振动响应','Color','black','fontsize',10);
rectangle('position',[5,8,1,0.1],'FaceColor','g');
text(0,8.1,'质量块3振动响应','Color','black','fontsize',10);
title('模态坐标系-解析解');%S=M^(-1/2)*P,x=S*y,此时求解为y坐标下响应,若求物理坐标响应,则需要乘S
t1 = linspace(0,30,1000);
C1=double(subs(C1,t,t1));
D1=double(subs(D1,t,t1));
G1=double(subs(G1,t,t1));
for k = 1:length(t1)
    %first line
    y= C1;
   addpoints(h1,t1(k),y(k));
   
   %second line
   y=D1;
   addpoints(h2,t1(k),y(k));
   
   %third line
   y=G1;
   addpoints(h3,t1(k),y(k));
   
   %update screen
   drawnow
end

x=S*[r1;r2;r3];
x1=x(1,1);
x2=x(2,1);
x3=x(3,1);
x(1)=simplify(x1);
x(2)=simplify(x2);
x(3)=simplify(x3);
C2=0;D2=0;G2=0;
subplot(1,2,2);
C2=eval(x(1));
D2=eval(x(2));
G2=eval(x(3));
C2=double(subs(C2,t,t1));
D2=double(subs(D2,t,t1));
G2=double(subs(G2,t,t1));
h1=animatedline('Color','b','LineStyle','-','LineWidth',2);
h2=animatedline('Color','r','LineStyle','-.','LineWidth',2);
h3=animatedline('Color','g','LineStyle','--','LineWidth',2);
axis([0,30,-15,10])
set(gca,'XTick',0:30/5:30); %设定x轴坐标刻度
set(gca,'YTick',-15:25/5:10); %设定y轴坐标刻度
xlabel('时间(T/s)')
rectangle('position',[5,9,1,0.1],'FaceColor','b');
text(0,9.1,'质量块1振动响应','Color','black','fontsize',10);
rectangle('position',[5,8.5,1,0.1],'FaceColor','r');
text(0,8.6,'质量块2振动响应','Color','black','fontsize',10);
rectangle('position',[5,8,1,0.1],'FaceColor','g');
text(0,8.1,'质量块3振动响应','Color','black','fontsize',10);
title('物理坐标系-解析解');
t1 = linspace(0,30,1000);
for k = 1:length(t1)
    %first line
    y= C2;
   addpoints(h1,t1(k),y(k));
   
   %second line
   y=D2;
   addpoints(h2,t1(k),y(k));
   
   %third line
   y=G2;
   addpoints(h3,t1(k),y(k));
   
   %update screen
   drawnow
end
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

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