关于增量谐波平衡法（IHB）的编程问题

如枫

有没有人编过增量谐波平衡法的Mathematica程序？我最近在编程序，一直都不行，想请教一下
下面是我编的Van der pol方程程序，但是好像不收敛，不知道是程序有问题还是初值选取不好

1. Unprotect["`*"];
   
2. ClearAll["`*"];
   
3. (*输入初始数据*)
   
4. nc = 1; ns = 1; n = nc + ns;
   
5. omega0 = 0.94;
   
6. epsilon = 1.0;
   
7. (*输入初始矩阵*)
   
8. cba = epsilon;          (*输入C拔矩阵*)
   
9. kba = 1;      (*输入K拔矩阵*)
   
10. 
    
11. a = Array[a0, n];      (*输入A的初始矩阵*)
    
12. a[[nc + 1]] = 2;
    
13. Do[a[[i]] = 0.241, {i, 1, nc}];
    
14. Do[a[[i]] = 0.126, {i, nc + 2, n}];
    
15. 
    
16. (*主程序*)
    
17. s = Array[cs0, {1, n}];
    
18. Do[cs0[1, i] = Cos[i*t], {i, nc}];
    
19. Do[cs0[1, i] = Sin[(i - nc)*t], {i, nc + 1, n}];
    
20. 
    
21. s1 = Dt[s, t];
    
22. s2 = Dt[s1, t];
    
23. st = Transpose[s];
    
24. m0 = st.s2;
    
25. m = Integrate[m0, {t, 0, 2*Pi}];
    
26. c0 = cba st.s1;
    
27. c = Integrate[c0, {t, 0, 2*Pi}];
    
28. k0 = kba st.s;
    
29. k = Integrate[k0, {t, 0, 2*Pi}];
    
30. 
    
31. r = Array[r0, {n, 1}];
    
32. Do[r0[i, 1] = 1, {i, n}];
    
33. maxr = 1;
    
34. nn = 1;
    
35. x = s.a;
    
36. c3ba = epsilon x^2;
    
37. x1 = Dt[x, t];
    
38. k2ba = epsilon x x1;
    
39. c30 = c3ba[[1]] st.s1;
    
40. k20 = k2ba[[1]] st.s;
    
41. c3 = Integrate[c30, {t, 0, 2*Pi}];
    
42. k2 = Integrate[k20, {t, 0, 2*Pi}];
    
43. kmc = omega0^2 m + omega0 (c3 - c) + k + 2 omega0 k2;
    
44. rmc = -(2 omega0 m + c3 - c).a;
    
45. r = -(omega0^2 m + omega0 (c3 - c) + k + k2).a;
    
46. While[maxr > 0.00001,
    
47. 
    
48. (*将kmc中第一列元素换为rmc*)
    
49. 
    
50. 
    
51. Do[kmc[[i, nc + 1]] = -rmc[[i]], {i, 1, n}];
    
52. 
    
53. 
    
54. 
    
55. (*定义方程*)
    
56. (*left=kmc.da;
    
57. right=r-domega*rmc;*)
    
58. 
    
59. 
    
60. (*If[Det[kmc]<0.1*10^(-10),
    
61. omega0=omega0+domega(*如果omega0增大不行，改成减小*);Continue[]];*)
    
62. da = LinearSolve[kmc, r, ZeroTest -> (PossibleZeroQ[Simplify[#]] &)];
    
63. (*Print["a0[1,1]="]
    
64. a0[1,1]
    
65. a0[1,1]+da[[1,1]]
    
66. omega0
    
67. omega0+domega*)
    
68. omega0 = omega0 + da[[nc + 1]];
    
69. Do[a[[i]] = a[[i]] + da[[i]], {i, 1, nc}];
    
70. Do[a[[i]] = a[[i]] + da[[i]], {i, nc + 2, n}];
    
71. (*omega0=omega0+domega;*)
    
72. (*a0[1,1]
    
73. omega0*)
    
74. x = s.a;
    
75. c3ba = epsilon x^2;
    
76. x1 = Dt[x, t];
    
77. k2ba = epsilon x x1;
    
78. c30 = c3ba[[1]] st.s1;
    
79. k20 = k2ba[[1]] st.s;
    
80. c3 = Integrate[c30, {t, 0, 2*Pi}];
    
81. k2 = Integrate[k20, {t, 0, 2*Pi}];
    
82. kmc = omega0^2 m + omega0 (c3 - c) + k + 2 omega0 k2;
    
83. rmc = -(2 omega0 m + c3 - c).a;
    
84. r = -(omega0^2 m + omega0 (c3 - c) + k + k2).a;
    
85. maxr = Max[Abs[Take[r, n]]];
    
86. nn = nn + 1;
    
87. If[nn > 10, Print["求解失败"]; Break[]];
    
88. ]
    
89. 
    
90. Print["x=" MatrixForm[s.a]];
    
91. Print["Maxr="];
    
92. maxr
    
93. Print["Det[kmc]="];
    
94. Det[kmc]
    
95. Print["da"];
    
96. da

复制代码

charminggu

楼主能不能给一个具体一点的增量谐波平衡法的定义呀？

如枫

楼上可以看看陈树辉的《强非线性振动系统的定量分析方法》，编程倒问题不大，关键是选初始解太难了，
我请教过陈树辉的学生，也说没什么好的办法！

xuruikl

兄弟，你说的陈树辉《强非线性振动系统的定量分析方法》是电子版的吗？

如枫

有纸质版的，没见到过电子版的

zhangwenjing

收敛怎么解决？

kangarooli

回复 如枫 的帖子

程序一般没什么问题，主要就是取值，最近我也在做，用的是多尺度法，也是总仿不出来，我觉得可能还是没取对值
   

liuwg14

回复 kangarooli 的帖子

请问,你问题解决了吗?如何选取初值,有什么技巧没有呀?

kangarooli

回复 liuwg14 的帖子

我的也没解决，也没什么技巧，我就是多取几个试的

liuwg14

谢谢,我发现初值太敏感了.稍微有点改动,就得不到解. 我求解的是Mathieu方程!

liuwg14

这个问题不解决,后期工作不好弄呀

如枫

大家讨论一下初值的问题啊！能不能改进一下饭算法，使程序对初值不要这么敏感。
感觉数值解对初值选取有点帮助！

雪缘

如枫 发表于 2010-9-26 20:52

                               
登录/注册后可看大图

大家讨论一下初值的问题啊！能不能改进一下饭算法，使程序对初值不要这么敏感。
感觉数值解对初值选取有点 ...

这个问题估计不好办，对初值是否敏感是有方程自身的特性决定的
如果改变了那还是原来的方程吗？

个人认为还是从寻找初值的范围入手更加合理点

如枫

自治系统好办一点，非自治系统的程序难写，总是不收敛！有好多文献用IHB法求非自治系统的解，不知道具体怎么操作的！

Vickyvictoria

如枫 发表于 2010-10-4 20:15

                               
登录/注册后可看大图

自治系统好办一点，非自治系统的程序难写，总是不收敛！有好多文献用IHB法求非自治系统的解，不知道具体怎么 ...

我这儿有一个matlab版的IHB法程序，网站发现的，下载了还没用过，贴出来大家看看是否好使

1. function y=func()
   
2. clc
   
3. clear
   
4. syms t a0 a1 a2 a3 a4 a5 a6 a7 a8 b1 b2 b3 b4 b5 b6 b7 b8;
   
5. 
   
6. epsilon=0.125;
   
7. alpha=1.0;
   
8. omega=2.0;
   
9. gamma=1.0;
   
10. beta=4.6;
    
11. 
    
12. x=a0+a1*cos(t)+a2*cos(2*t)+b1*sin(t)+b2*sin(2*t);
    
13. %x=a0+a1*cos(t)+a2*cos(2*t)+a3*cos(3*t)+a4*cos(4*t)+a5*cos(5*t)+a6*cos(6*t)+a7*cos(7*t)+a8*cos(8*t)+b1*sin(t)+b2*sin(2*t)+b3*sin(3*t)+b4*sin(4*t)+b5*sin(5*t)+b6*sin(6*t)+b7*sin(7*t)+b8*sin(8*t);
    
14. 
    
15. % f=diff(diff(x,t),t)+2*epsilon*diff(x,t)-(alpha+beta*sin(omega*t))*x-gamma*x^3;
    
16. %f=omega^2*diff(diff(x,t),t)+2*omega*epsilon*diff(x,t)-(alpha+beta*sin(t))*x-gamma*x^3
    
17. d_x=diff(x,'t');
    
18. dd_x=diff(d_x,'t');
    
19. 
    
20. f=omega^2*dd_x+2*epsilon*omega*d_x-(alpha+beta*sin(t))*x+gamma*x^3
    
21. 
    
22. f_1=simplify(int(f,t,0,2*pi)/(2*pi));
    
23. f_cost=simplify(int(f*cos(t),t,0,2*pi)/pi);
    
24. f_cos2t=simplify(int(f*cos(2*t),t,0,2*pi)/pi);
    
25. f_sint=simplify(int(f*sin(t),t,0,2*pi)/pi);
    
26. f_sin2t=simplify(int(f*sin(2*t),t,0,2*pi)/pi);
    
27. 
    
28. % fun=[f_1;f_cost;f_cos2t;f_sint;f_sin2t]
    
29. % x0(5)=0
    
30. % Newtons(fun,x0)

复制代码