增量谐波平衡法 IHB FDB

nedusts

　　大谢本程序的编写者和分享者!!!!!!
　　与程序相关的有用的参考文献：《应用于范德波方程的增量谐波平衡法》唐南中山大学研究生专刊自然科学版

1. 　　clear all
   
2. 
   
3. 　　clc
   
4. 
   
5. 　　tic
   
6. 
   
7. 　　%%%定义各参数
   
8. 
   
9. 　　syms t
   
10. 
    
11. 　　w0=3;
    
12. 
    
13. 　　epsR=0.001;
    
14. 
    
15. 　　m=[1 0;0 1];
    
16. 
    
17. 　　epsilon=0.24;
    
18. 
    
19. 　　r=0.4;
    
20. 
    
21. 　　delta=0.56;
    
22. 
    
23. 　　k=[1+epsilon*r -epsilon*r;-epsilon*r 1+epsilon*delta];
    
24. 
    
25. 　　Cs=[cos(t) cos(3*t) sin(t) sin(3*t)];
    
26. 
    
27. 　　Cs1=diff(Cs,t,1);
    
28. 
    
29. 　　S=[cos(t) cos(3*t) sin(t) sin(3*t) 0 0 0 0;0 0 0 0 cos(t) cos(3*t) sin(t) sin(3*t)];
    
30. 
    
31. 　　A1=[1 1 1 1]';
    
32. 
    
33. 　　A2=[1 1 1 1]';
    
34. 
    
35. 　　A0=[A1;A2];
    
36. 
    
37. 　　T1=[eye(4,4) zeros(4,4)];
    
38. 
    
39. 　　T2=[zeros(4,4) eye(4,4)];
    
40. 
    
41. 　　S2=diff(S,t,2);
    
42. 
    
43. 　　fm=inline(S'*m*S2);
    
44. 
    
45. 　　M=quadv(fm,0,2*pi);
    
46. 
    
47. 　　fk=inline(S'*k*S);
    
48. 
    
49. 　　K=quadv(fk,0,2*pi);
    
50. 
    
51. 　　S1=diff(S,t,1);
    
52. 
    
53. 　　c=[1 0;0 1];
    
54. 
    
55. 　　fc=inline(S'*c*S1);
    
56. 
    
57. 　　C=quadv(fc,0,2*pi);
    
58. 
    
59. 　　c3=diag(S*A0).^2;
    
60. 
    
61. 　　fc3=inline(S'*c3*S1);
    
62. 
    
63. 　　C3=quadv(fc3,0,2*pi);
    
64. 
    
65. 　　k2=diag(S*A0).*diag(S1*A0)
    
66. 
    
67. 　　fk2=inline(S'*k2*S);
    
68. 
    
69. 　　K2=quadv(fk2,0,2*pi);
    
70. 
    
71. 　　%%%%%带入推导出的公式
    
72. 
    
73. 　　Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;
    
74. 
    
75. 　　R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;
    
76. 
    
77. 　　Rmc=-(2*w0*M+epsilon*(C3-C))*A0;
    
78. 
    
79. 　　% AA首元素已知a1=0.0,求ww
    
80. 
    
81. 　　a1=0.0;
    
82. 
    
83. 　　%%%变换矩阵，使ww变量代替a1
    
84. 
    
85. 　　Kmc11=-Rmc(:,1);
    
86. 
    
87. 　　Kmcr=[Kmc11 Kmc(:,2:size(Kmc,2))];
    
88. 
    
89. 　　%求未知变量
    
90. 
    
91. 　　AA=inv(Kmcr)*R; % drtA1(1)
    
92. 
    
93. 　　ww=AA(1); % drtW(1)
    
94. 
    
95. 　　%%赋予新变量新值
    
96. 
    
97. 　　A01=A0+[a1;AA(2:length(A0),1)]; % A(1)+drtA(1)
    
98. 
    
99. 　　% Aw0=AA+A00; % A1(0)+drtA1(1)=A1(1)
    
100. 
     
101. 　　w01=w0+ww; % W+drtW(1)
     
102. 
     
103. 　　n=1;
     
104. 
     
105. 　　tol=1;
     
106. 
     
107. 　　%
     
108. 
     
109. 　　while tol>epsR
     
110. 
     
111. 　　A0=A01;
     
112. 
     
113. 　　%
     
114. 
     
115. 　　w0=w01
     
116. 
     
117. 　　%
     
118. 
     
119. 　　c3=diag(S*A0).^2;
     
120. 
     
121. 　　fc3=inline(S'*c3*S1);
     
122. 
     
123. 　　C3=quadv(fc3,0,2*pi);
     
124. 
     
125. 　　k2=diag(S*A0).*diag(S1*A0)
     
126. 
     
127. 　　fk2=inline(S'*k2*S);
     
128. 
     
129. 　　K2=quadv(fk2,0,2*pi);
     
130. 
     
131. 　　%%%%%带入推导出的公式
     
132. 
     
133. 　　Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;
     
134. 
     
135. 　　R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;
     
136. 
     
137. 　　Rmc=-(2*w0*M+epsilon*(C3-C))*A0;
     
138. 
     
139. 　　%%%%%
     
140. 
     
141. 　　tol=norm(R)
     
142. 
     
143. 　　if(n>1000)
     
144. 
     
145. 　　disp('迭代步数太多，可能不收敛')
     
146. 
     
147. 　　return;
     
148. 
     
149. 　　end
     
150. 
     
151. 　　Kmc11=-Rmc(:,1);
     
152. 
     
153. 　　Kmcr=[Kmc11 Kmc(:,2:size(Kmc,2))];
     
154. 
     
155. 　　%
     
156. 
     
157. 　　AA=inv(Kmcr)*R;
     
158. 
     
159. 　　ww=AA(1);
     
160. 
     
161. 　　% A00=[w0;A0(2:6,1)];
     
162. 
     
163. 　　A01=A0+[a1;AA(2:length(A0),1)];
     
164. 
     
165. 　　w01=w0+ww;
     
166. 
     
167. 　　n=n+1;
     
168. 
     
169. 　　end
     
170. 
     
171. 　　X0=S*A0;
     
172. 
     
173. 　　dX0=S1*A0;
     
174. 
     
175. 　　%%绘范德波图
     
176. 
     
177. 　　tt=0:.1:10;
     
178. 
     
179. 　　xo1=subs(X0(1),tt);
     
180. 
     
181. 　　xo2=subs(X0(2),tt);
     
182. 
     
183. 　　dxo1=subs(dX0(1),tt);
     
184. 
     
185. 　　dxo2=subs(dX0(2),tt);
     
186. 
     
187. 　　figure(1)
     
188. 
     
189. 　　plot(xo1,dxo1,'b','linewidth',2)
     
190. 
     
191. 　　hold on
     
192. 
     
193. 　　plot(xo2,dxo2,'b','linewidth',2)
     
194. 
     
195. 　　axis([-3 3 -3 3])
     
196. 
     
197. 　　title('范德波极限环')
     
198. 
     
199. 　　xlabel('x0')
     
200. 
     
201. 　　ylabel('dx0')
     
202. 
     
203. 　　toc
     

复制代码

转自：http://blog.sina.com.cn/s/blog_786ce14d0102v1pa.html

lihaitao123

请问一下，A1=[1 1 1 1]'; 　　A2=[1 1 1 1]'; 是表示初值么？ 可以改变么？

花花三公子

lihaitao123 发表于 2017-5-20 21:32
请问一下，A1=[1 1 1 1]'; 　　A2=[1 1 1 1]'; 是表示初值么？ 可以改变么？

是初值，可以修改

花花三公子

这个程序太早了，用的inline和quadv函数，这些函数再新版matlab里会提示：未来版本中会移除