设为首页收藏本站联系我们
切换到宽版

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛»论坛 › 基础理论 › 数学与非线性科学 › 关于matlab求解duffing
返回列表 发新帖
查看: 2494|回复: 6

[分形与混沌] 关于matlab求解duffing

[复制链接]
不是一样的
发表于 2017-3-9 15:39 | 显示全部楼层 |阅读模式

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要 登录 才可以下载或查看，没有账号？我要加入

x
小白一枚  怎么用matlab求解duffing方程   求程序  有解释更好了  谢谢大家
回复
分享到：

使用道具 举报

Agoni
发表于 2017-3-10 09:16 | 显示全部楼层
http://forum.vibunion.com/thread-151083-1-1.html
参考这个
回复 支持 反对

使用道具 举报

不是一样的
 楼主| 发表于 2017-3-10 15:45 | 显示全部楼层
看不懂啊   有详细的步骤吗   有时间的话希望详细解答
回复 支持 反对

使用道具 举报

kuzb
发表于 2017-3-15 11:36 | 显示全部楼层
“求解”指的是什么
回复 支持 反对

使用道具 举报

truleeee
发表于 2017-3-15 14:08 | 显示全部楼层
  1. function dx=Duffing(t,x,flag,f); global f; %Duffing方程
  2. a=0.25;
  3. b=-1;
  4. c=1;
  5. w=1;
  6. dx=[x(2);-a*x(2)-b*x(1)-c*x(1)^3+f*cos(x(3));w];
  7. %fftplot幅值谱
  8. function [f,y]=fftplot(x,step)
  9. [r,c] = size(x);
  10. if r == 1
  11.         x = x.';   % make it a column
  12. end;
  13. [n,cc] = size(x);
  14. m = 2^nextpow2(n);
  15. y = fft(x,m);
  16. y = y(1:n,:);
  17. yy=abs(y)*step;
  18. fid1=fopen('fft1.txt','w');
  19. fid2=fopen('aft1.txt','w');
  20. for kk=1:m/2+1;
  21.    f(kk)=(kk-1)/m/step;
  22.    fprintf(fid1,'%8.4f\n',f(kk));
  23.    fprintf(fid2,'%8.4f\n',y(kk));
  24. end
  25. fclose(fid1);
  26. fclose(fid2);
  27. y=yy(1:m/2+1);
  28. %plot(f,y);
  29. %xlabel('Frequency (Hz)');
  30. %ylabel('Fourier Amplitude');
  31. %title('Fouier Transform of the Data')
  32. %求解、绘图：时域波形、相平面、幅值谱、Poincare图、分岔图
  33. clc
  34. global f;
  35. f=0.1
  36. tspan=[0:0.01*2*pi:2000];
  37. options=odeset('RelTol',1e-5,'AbsTol',1e-6);
  38. [t,x]=ode45('Duffing',tspan,[1;1.5;0],options);
  39. x_data=x(:,1);
  40. x_dot=x(:,2);
  41. %时域波形
  42. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  43. hold on;
  44. plot(tspan,x_data,'LineWidth',1.5);
  45. xlim([100 200])
  46. %xlim([200 300])
  47. xlabel('\itt','Fontsize',20,'Fontname','Times new roman');%x轴标注
  48. ylabel('\itx','Fontsize',20,'Fontname','Times new roman','Rotation',0);%y轴标注
  49. grid on
  50. box on
  51. %相平面图
  52. figure(2)
  53. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  54. hold on;
  55. plot(x_data(1000:3000),x_dot(1000:3000),'LineWidth',1.5);
  56. %plot(x_data(10:30000),x_dot(10:30000),'LineWidth',1.5);
  57. xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
  58. ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
  59. grid on
  60. box on
  61. %幅值谱
  62. [freq,amp]=fftplot(x_data,0.01*2*pi);
  63. figure(3)
  64. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  65. hold on;
  66. xlim([0 1])
  67. ylim([0 200])
  68. plot(freq(1:5000),amp(1:5000),'LineWidth',1.5);
  69. xlabel('/Hz','Fontsize',20,'Fontname','Times new roman');%x轴标注
  70. ylabel('mu','Fontsize',20,'Fontname','Times new roman');%y轴标注
  71. grid on
  72. box on
  73. %poincare图
  74. x(:,2)=mod(x(:,2),2*pi)-pi;
  75. phi0=pi*1;
  76. for k=1:round(max(x(:,3))/2/pi);
  77.     d=x(:,3)-(k-1)*2*pi-phi0;
  78.     [P,K]=sort(abs(d));
  79.     x1l=x(K(1),1);
  80.     x1r=x(K(2),1);
  81.     x2l=x(K(1),2);
  82.     x2r=x(K(2),2);
  83.     x3l=x(K(1),3);
  84.     x3r=x(K(2),3);
  85.     if abs(P(1))+abs(P(2))<3e-16;
  86.         X1(k)=x1l;
  87.         X2(k)=x2l;
  88.     else
  89.         Q=polyfit([x3l,x3r],[x1l,x1r],1);
  90.         X1(k)=polyval(Q,(k-1)*2*pi-phi0);
  91.          Q=polyfit([x3l,x3r],[x2l,x2r],1);
  92.         X2(k)=polyval(Q,(k-1)*2*pi-phi0);
  93.     end
  94. end
  95.     figure(4)
  96.    axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  97.    plot(X1,X2,'.');
  98.   xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
  99. ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
  100. %分岔图
  101. global f
  102. figure(5)
  103. hold on;
  104. box on;
  105. xlim([0,1]);
  106. ylim([-2,2]);
  107. N=80;
  108. Tn=zeros(1,N);
  109. options=odeset('RelTol',1e-9);
  110. for f=0.05:.05:1;
  111. x0=[0,1.5,0];
  112. for n=1:60;
  113. [t,x]=ode45('Duffing',[0,2*pi],x0,options);
  114. x0=x(end,:);
  115. end
  116. for n=1:N;
  117. [t,x]=ode45('Duffing',[0,pi*2],x0,options);
  118. x0=x(end,:);
  119. xd=x(:,1);
  120. Tn(n)=max(xd);
  121. end
  122. f
  123. plot(f,Tn,'b.','markersize',1);
  124. pause(0.0001);
  125. end
  126. xlabel('\itf','Fontsize',16,'Fontname','Times new roman');%x轴标注
  127. ylabel('\itx','Fontsize',16,'Fontname','Times new roman','Rotation',0);%y轴标注
复制代码


回复 支持 反对

使用道具 举报

truleeee
发表于 2017-3-15 14:08 | 显示全部楼层
  1. function dx=Duffing(t,x,flag,f); global f; %Duffing方程
  2. a=0.25;
  3. b=-1;
  4. c=1;
  5. w=1;
  6. dx=[x(2);-a*x(2)-b*x(1)-c*x(1)^3+f*cos(x(3));w];
  7. %fftplot幅值谱
  8. function [f,y]=fftplot(x,step)
  9. [r,c] = size(x);
  10. if r == 1
  11.         x = x.';   % make it a column
  12. end;
  13. [n,cc] = size(x);
  14. m = 2^nextpow2(n);
  15. y = fft(x,m);
  16. y = y(1:n,:);
  17. yy=abs(y)*step;
  18. fid1=fopen('fft1.txt','w');
  19. fid2=fopen('aft1.txt','w');
  20. for kk=1:m/2+1;
  21.    f(kk)=(kk-1)/m/step;
  22.    fprintf(fid1,'%8.4f\n',f(kk));
  23.    fprintf(fid2,'%8.4f\n',y(kk));
  24. end
  25. fclose(fid1);
  26. fclose(fid2);
  27. y=yy(1:m/2+1);
  28. %plot(f,y);
  29. %xlabel('Frequency (Hz)');
  30. %ylabel('Fourier Amplitude');
  31. %title('Fouier Transform of the Data')
  32. %求解、绘图：时域波形、相平面、幅值谱、Poincare图、分岔图
  33. clc
  34. global f;
  35. f=0.1
  36. tspan=[0:0.01*2*pi:2000];
  37. options=odeset('RelTol',1e-5,'AbsTol',1e-6);
  38. [t,x]=ode45('Duffing',tspan,[1;1.5;0],options);
  39. x_data=x(:,1);
  40. x_dot=x(:,2);
  41. %时域波形
  42. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  43. hold on;
  44. plot(tspan,x_data,'LineWidth',1.5);
  45. xlim([100 200])
  46. %xlim([200 300])
  47. xlabel('\itt','Fontsize',20,'Fontname','Times new roman');%x轴标注
  48. ylabel('\itx','Fontsize',20,'Fontname','Times new roman','Rotation',0);%y轴标注
  49. grid on
  50. box on
  51. %相平面图
  52. figure(2)
  53. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  54. hold on;
  55. plot(x_data(1000:3000),x_dot(1000:3000),'LineWidth',1.5);
  56. %plot(x_data(10:30000),x_dot(10:30000),'LineWidth',1.5);
  57. xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
  58. ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
  59. grid on
  60. box on
  61. %幅值谱
  62. [freq,amp]=fftplot(x_data,0.01*2*pi);
  63. figure(3)
  64. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  65. hold on;
  66. xlim([0 1])
  67. ylim([0 200])
  68. plot(freq(1:5000),amp(1:5000),'LineWidth',1.5);
  69. xlabel('/Hz','Fontsize',20,'Fontname','Times new roman');%x轴标注
  70. ylabel('mu','Fontsize',20,'Fontname','Times new roman');%y轴标注
  71. grid on
  72. box on
  73. %poincare图
  74. x(:,2)=mod(x(:,2),2*pi)-pi;
  75. phi0=pi*1;
  76. for k=1:round(max(x(:,3))/2/pi);
  77.     d=x(:,3)-(k-1)*2*pi-phi0;
  78.     [P,K]=sort(abs(d));
  79.     x1l=x(K(1),1);
  80.     x1r=x(K(2),1);
  81.     x2l=x(K(1),2);
  82.     x2r=x(K(2),2);
  83.     x3l=x(K(1),3);
  84.     x3r=x(K(2),3);
  85.     if abs(P(1))+abs(P(2))<3e-16;
  86.         X1(k)=x1l;
  87.         X2(k)=x2l;
  88.     else
  89.         Q=polyfit([x3l,x3r],[x1l,x1r],1);
  90.         X1(k)=polyval(Q,(k-1)*2*pi-phi0);
  91.          Q=polyfit([x3l,x3r],[x2l,x2r],1);
  92.         X2(k)=polyval(Q,(k-1)*2*pi-phi0);
  93.     end
  94. end
  95.     figure(4)
  96.    axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  97.    plot(X1,X2,'.');
  98.   xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
  99. ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
  100. %分岔图
  101. global f
  102. figure(5)
  103. hold on;
  104. box on;
  105. xlim([0,1]);
  106. ylim([-2,2]);
  107. N=80;
  108. Tn=zeros(1,N);
  109. options=odeset('RelTol',1e-9);
  110. for f=0.05:.05:1;
  111. x0=[0,1.5,0];
  112. for n=1:60;
  113. [t,x]=ode45('Duffing',[0,2*pi],x0,options);
  114. x0=x(end,:);
  115. end
  116. for n=1:N;
  117. [t,x]=ode45('Duffing',[0,pi*2],x0,options);
  118. x0=x(end,:);
  119. xd=x(:,1);
  120. Tn(n)=max(xd);
  121. end
  122. f
  123. plot(f,Tn,'b.','markersize',1);
  124. pause(0.0001);
  125. end
  126. xlabel('\itf','Fontsize',16,'Fontname','Times new roman');%x轴标注
  127. ylabel('\itx','Fontsize',16,'Fontname','Times new roman','Rotation',0);%y轴标注
复制代码


回复 支持 反对

使用道具 举报

truleeee
发表于 2017-3-15 14:09 | 显示全部楼层
function dx=Duffing(t,x,flag,f); global f; %Duffing方程
a=0.25;
b=-1;
c=1;
w=1;
dx=[x(2);-a*x(2)-b*x(1)-c*x(1)^3+f*cos(x(3));w];
%fftplot幅值谱
function [f,y]=fftplot(x,step)
[r,c] = size(x);
if r == 1
        x = x.';   % make it a column
end;
[n,cc] = size(x);
m = 2^nextpow2(n);
y = fft(x,m);
y = y(1:n,:);
yy=abs(y)*step;
fid1=fopen('fft1.txt','w');
fid2=fopen('aft1.txt','w');
for kk=1:m/2+1;
   f(kk)=(kk-1)/m/step;
   fprintf(fid1,'%8.4f\n',f(kk));
   fprintf(fid2,'%8.4f\n',y(kk));
end
fclose(fid1);
fclose(fid2);
y=yy(1:m/2+1);
%plot(f,y);
%xlabel('Frequency (Hz)');
%ylabel('Fourier Amplitude');
%title('Fouier Transform of the Data')
%求解、绘图：时域波形、相平面、幅值谱、Poincare图、分岔图
clc
global f;
f=0.1
tspan=[0:0.01*2*pi:2000];
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
[t,x]=ode45('Duffing',tspan,[1;1.5;0],options);
x_data=x(:,1);
x_dot=x(:,2);
%时域波形
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
hold on;
plot(tspan,x_data,'LineWidth',1.5);
xlim([100 200])
%xlim([200 300])
xlabel('\itt','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('\itx','Fontsize',20,'Fontname','Times new roman','Rotation',0);%y轴标注
grid on
box on
%相平面图
figure(2)
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
hold on;
plot(x_data(1000:3000),x_dot(1000:3000),'LineWidth',1.5);
%plot(x_data(10:30000),x_dot(10:30000),'LineWidth',1.5);
xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
grid on
box on
%幅值谱
[freq,amp]=fftplot(x_data,0.01*2*pi);
figure(3)
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
hold on;
xlim([0 1])
ylim([0 200])
plot(freq(1:5000),amp(1:5000),'LineWidth',1.5);
xlabel('/Hz','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('mu','Fontsize',20,'Fontname','Times new roman');%y轴标注
grid on
box on
%poincare图
x(:,2)=mod(x(:,2),2*pi)-pi;
phi0=pi*1;
for k=1:round(max(x(:,3))/2/pi);
    d=x(:,3)-(k-1)*2*pi-phi0;
    [P,K]=sort(abs(d));
    x1l=x(K(1),1);
    x1r=x(K(2),1);
    x2l=x(K(1),2);
    x2r=x(K(2),2);
    x3l=x(K(1),3);
    x3r=x(K(2),3);
    if abs(P(1))+abs(P(2))<3e-16;
        X1(k)=x1l;
        X2(k)=x2l;
    else
        Q=polyfit([x3l,x3r],[x1l,x1r],1);
        X1(k)=polyval(Q,(k-1)*2*pi-phi0);
         Q=polyfit([x3l,x3r],[x2l,x2r],1);
        X2(k)=polyval(Q,(k-1)*2*pi-phi0);
    end
end
    figure(4)
   axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
   plot(X1,X2,'.');
  xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
%分岔图
global f
figure(5)
hold on;
box on;
xlim([0,1]);
ylim([-2,2]);
N=80;
Tn=zeros(1,N);
options=odeset('RelTol',1e-9);
for f=0.05:.05:1;
x0=[0,1.5,0];
for n=1:60;
[t,x]=ode45('Duffing',[0,2*pi],x0,options);
x0=x(end,:);
end
for n=1:N;
[t,x]=ode45('Duffing',[0,pi*2],x0,options);
x0=x(end,:);
xd=x(:,1);
Tn(n)=max(xd);
end
f
plot(f,Tn,'b.','markersize',1);
pause(0.0001);
end
xlabel('\itf','Fontsize',16,'Fontname','Times new roman');%x轴标注
ylabel('\itx','Fontsize',16,'Fontname','Times new roman','Rotation',0);%y轴标注
回复 支持 反对

使用道具 举报

返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-19 22:19 , Processed in 0.064658 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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