[求助]求随机共振源程序

xcsyb

谢谢各位！！！

风花雪月

随机共振的基本含义是指一个非线性双稳或多稳系统，当在小周期信号驱动下不足使系统的输出在两个稳态之间跳跃，只能在一个稳态上。而在弱噪声和小周期调制信号共同作用下，随着输入噪声强度的增加，系统输出功率谱中调制信号的频率处出现一个峰值．当噪声强度达到某一适当值时，系统在两定态之间进行跃迁，输出信号的峰值达到最大——“共振”。以后随着噪声强度的增强，输出信号的峰值逐渐减小。总结产生随机共振现象的条件，大致可以归纳为以下三个要素： <BR>（１）非线性系统要存在两个（或多个）稳态。 <BR>（２）输入信号，并且输入信号还不足引起系统在不同稳态之间跃迁。 <BR>（３）噪声。它既可以是系统外部的随机驱动力，也可以来自系统内部相互作用引起的涨落。 <BR>

风花雪月

下面是matlab写的一个随机共振的例子，你可以参考一下

1. %programshijin.m
   
2. x=-2:0.1:2;区间范围
   
3. b=1;势阱参数,决定势阱的深度与宽度
   
4. c=1;
   
5. 
   
6. figure(1)
   
7. subplot(2,1,1);分区作图
   
8. y=-b.*x.^2/2+c.*x.^4/4;势函数
   
9. plot(x,y);作二维图形
   
10. subplot(2,1,2);
    
11. b=3:0.1:7;给定b,c的变化区间
    
12. c=2:0.1:6;
    
13. h=moviein(40);影片动画
    
14. fori=1:40;
    
15. y=-b(i).*x.^2/2+c(i).*x.^4/4;
    
16. plot(x,y);
    
17. h(:,i)=getframe;
    
18. end
    
19. title('势函数随b,c的变化情况')坐标标注
    
20. subplot(2,1,1)
    
21. title('无信号及噪声输入时系统的势函数')
    
22. xlabel('x');
    
23. ylabel('U(x)');
    
24. subplot(2,1,2)
    
25. xlabel('x');
    
26. ylabel('U(x)');
    
27. movie(h,1);
    
28. figure(2)
    
29. b=1;给定b,c加入周期调制信号
    
30. c=1;
    
31. w=pi;
    
32. z=0;
    
33. subplot(2,1,1)信号较小时
    
34. A=2;
    
35. m=moviein(20);影片动画
    
36. [t,yy]=ode23('xinhao1',[0:0.1:10],[1,b-c+A],[],b,z,c,A,w);ode命令解微分方程
    
37. n=length(t);
    
38. fori=1:fix(n/2)
    
39. axis([-22-510]);
    
40. y=-b.*x.^2/2+c.*x.^4/4-A.*x.*sin(w*t(i));
    
41. plot(x,y);
    
42. holdon
    
43. uu=yy(i,2)*yy(i,2)/2-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
    
44. uu=-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
    
45. plot(yy(i,1),uu,'.r','markersize',10)
    
46. holdoff
    
47. m(:,i)=getframe;
    
48. end
    
49. title('加入周期信号较小，不足以翻转')
    
50. xlabel('x');
    
51. ylabel('U(x)');
    
52. v=moviein(20);
    
53. subplot(2,1,2)周期信号较大时
    
54. A=3;
    
55. [t,yy]=ode23('xinhao1',[0:0.1:10],[1,b-c+A],[],b,z,c,A,w);
    
56. n=length(t);
    
57. fori=1:fix(n/2)
    
58. y=-b.*x.^2/2+c.*x.^4/4-A.*x.*sin(w*t(i));
    
59. plot(x,y);
    
60. holdon
    
61. uu=yy(i,2)*yy(i,2)/2-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
    
62. uu=-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
    
63. plot(yy(i,1),uu,'.r','markersize',10)
    
64. holdoff
    
65. v(:,i)=getframe;
    
66. end
    
67. title('加入周期调制信号较大，可以翻转')
    
68. holdon
    
69. xlabel('x');
    
70. ylabel('U(x)');
    
71. axis([-22-510])
    
72. holdon
    
73. movie(v)
    
74. t=0:0.1:10;
    
75. n=length(t);
    
76. figure(3)加入信号同时加入随即噪声
    
77. A=2;
    
78. b=1;
    
79. m=moviein(20);
    
80. fori=1:fix(n/2)
    
81. r=2*rand(1)随即量
    
82. y=-b.*x.^2/2+c.*x.^4/4-A.*x.*sin(w*t(i))-r.*x;
    
83. plot(x,y);
    
84. holdon
    
85. [t,s,r]=oula('f',0,5,rand(1),0.1,r);调用欧拉函数解方程
    
86. u=-b.*s(i).^2/2+c.*s(i).^4/4-A.*s.*sin(w*t(i))-r*s(i);
    
87. plot(s(i),u(i),'.r','markersize',10);
    
88. holdoff
    
89. m(:,i)=getframe;
    
90. end
    
91. title('加入随机噪声的情况,参数选取同figure(2)中未翻转时');
    
92. xlabel('x');
    
93. ylabel('U(x)');
    
94. movie(m,1)
    
95. figure(4)波形图及频谱分析
    
96. x=-2:0.1:2;
    
97. w=50*pi;
    
98. b=1;
    
99. c=1
    
100. A=1;
     
101. t=0:0.01:10;
     
102. n=length(t);
     
103. subplot(2,1,1);
     
104. fori=1:n
     
105. y(i)=b.*x(1)-c.*x(1).^3+A.*cos(w*t(i))+0.5*rand(1);
     
106. end
     
107. plot(t,y)
     
108. title('信号的波形图');
     
109. xlabel('t');
     
110. ylabel('U(x)');
     
111. subplot(2,1,2)
     
112. z=fft(y);
     
113. z(1)=[];
     
114. n=length(z);
     
115. p=10^(-5).*abs(z(1:n/2)).^2
     
116. m=1/2;
     
117. f=100*(1:n/2)/(n/2)*m;
     
118. plot(f,p)
     
119. title('信号的频谱图,w=50*pi,f=25');
     
120. 
     
121. 
     
122. %oula.m子程序:欧拉法解方程
     
123. function[x,y,r]=oula(fun,x0,xf,y0,h,r)
     
124. n=fix((xf-x0)/h);
     
125. x(1)=x0;
     
126. y(1)=y0;
     
127. x(n)=0;
     
128. y(n)=0;
     
129. fori=1:(n-1)
     
130. x(i+1)=x0+i*h;
     
131. y1=y(i)+h*feval(fun,x(i),y(i),r);调用方程
     
132. y2=y(i)+h*feval(fun,x(i+1),y1,r);
     
133. y(i+1)=(y1+y2)/2;
     
134. end
     
135. 
     
136. 
     
137. %f.m方程文件
     
138. %参数b,c,A,w为1,故省略。
     
139. functionf=f(t,x,r)
     
140. f=x-x.^3+cos(t)-r;
     
141. 
     
142. 
     
143. %odefilexinhao1.mode文件
     
144. functionydot=xinhao1(t,x,flag,b,z,c,A,w)
     
145. ydot=[b*x(1)-c*x(1)^3+A*cos(w*t)+z;...
     
146. b*x(2)-3*c*x(1)*x(1)*x(2)-A*w*sin(w*t)];

复制代码

[ 本帖最后由 风花雪月 于 2008-4-30 10:30 编辑 ]

xcsyb

谢谢！！！

nanlaibeiqu

为什么我运行不了楼主的程序呢，提示
??? Error: File: sr1.m Line: 123 Column: 1
Function definitions are not permitted at the prompt or in scripts.

风花雪月

原帖由 nanlaibeiqu 于 2008-10-15 17:32 发表

                               
登录/注册后可看大图

为什么我运行不了楼主的程序呢，提示
??? Error: File: sr1.m Line: 123 Column: 1
Function definitions are not permitted at the prompt or in scripts.

都不知道你看过我上面的帖子没有，自己仔细看看吧，而不是盲目的粘贴运行

lihaitao123

风花雪月 发表于 2008-10-28 17:18

                               
登录/注册后可看大图

都不知道你看过我上面的帖子没有，自己仔细看看吧，而不是盲目的粘贴运行

谢谢风花雪月，我刚接触随机振动，太有帮助了。。

lordkibun

风花雪月 发表于 2006-6-25 09:19

                               
登录/注册后可看大图

下面是matlab写的一个随机共振的例子，你可以参考一下

[ 本帖最后由 风花雪月 于 2008-4-30 10:30 编辑 ]

您好。我也是在看这段程序。想问下，您用的是哪个版本的matlab？