声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 13741|回复: 25

[近似分析] 关于增量谐波平衡法(IHB)的编程问题

  [复制链接]
发表于 2009-2-20 08:36 | 显示全部楼层 |阅读模式

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

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

x
有没有人编过增量谐波平衡法的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. a = Array[a0, n];      (*输入A的初始矩阵*)
  11. a[[nc + 1]] = 2;
  12. Do[a[[i]] = 0.241, {i, 1, nc}];
  13. Do[a[[i]] = 0.126, {i, nc + 2, n}];

  14. (*主程序*)
  15. s = Array[cs0, {1, n}];
  16. Do[cs0[1, i] = Cos[i*t], {i, nc}];
  17. Do[cs0[1, i] = Sin[(i - nc)*t], {i, nc + 1, n}];

  18. s1 = Dt[s, t];
  19. s2 = Dt[s1, t];
  20. st = Transpose[s];
  21. m0 = st.s2;
  22. m = Integrate[m0, {t, 0, 2*Pi}];
  23. c0 = cba st.s1;
  24. c = Integrate[c0, {t, 0, 2*Pi}];
  25. k0 = kba st.s;
  26. k = Integrate[k0, {t, 0, 2*Pi}];

  27. r = Array[r0, {n, 1}];
  28. Do[r0[i, 1] = 1, {i, n}];
  29. maxr = 1;
  30. nn = 1;
  31. x = s.a;
  32. c3ba = epsilon x^2;
  33. x1 = Dt[x, t];
  34. k2ba = epsilon x x1;
  35. c30 = c3ba[[1]] st.s1;
  36. k20 = k2ba[[1]] st.s;
  37. c3 = Integrate[c30, {t, 0, 2*Pi}];
  38. k2 = Integrate[k20, {t, 0, 2*Pi}];
  39. kmc = omega0^2 m + omega0 (c3 - c) + k + 2 omega0 k2;
  40. rmc = -(2 omega0 m + c3 - c).a;
  41. r = -(omega0^2 m + omega0 (c3 - c) + k + k2).a;
  42. While[maxr > 0.00001,

  43. (*将kmc中第一列元素换为rmc*)


  44. Do[kmc[[i, nc + 1]] = -rmc[[i]], {i, 1, n}];



  45. (*定义方程*)
  46. (*left=kmc.da;
  47. right=r-domega*rmc;*)


  48. (*If[Det[kmc]<0.1*10^(-10),
  49. omega0=omega0+domega(*如果omega0增大不行,改成减小*);Continue[]];*)
  50. da = LinearSolve[kmc, r, ZeroTest -> (PossibleZeroQ[Simplify[#]] &)];
  51. (*Print["a0[1,1]="]
  52. a0[1,1]
  53. a0[1,1]+da[[1,1]]
  54. omega0
  55. omega0+domega*)
  56. omega0 = omega0 + da[[nc + 1]];
  57. Do[a[[i]] = a[[i]] + da[[i]], {i, 1, nc}];
  58. Do[a[[i]] = a[[i]] + da[[i]], {i, nc + 2, n}];
  59. (*omega0=omega0+domega;*)
  60. (*a0[1,1]
  61. omega0*)
  62. x = s.a;
  63. c3ba = epsilon x^2;
  64. x1 = Dt[x, t];
  65. k2ba = epsilon x x1;
  66. c30 = c3ba[[1]] st.s1;
  67. k20 = k2ba[[1]] st.s;
  68. c3 = Integrate[c30, {t, 0, 2*Pi}];
  69. k2 = Integrate[k20, {t, 0, 2*Pi}];
  70. kmc = omega0^2 m + omega0 (c3 - c) + k + 2 omega0 k2;
  71. rmc = -(2 omega0 m + c3 - c).a;
  72. r = -(omega0^2 m + omega0 (c3 - c) + k + k2).a;
  73. maxr = Max[Abs[Take[r, n]]];
  74. nn = nn + 1;
  75. If[nn > 10, Print["求解失败"]; Break[]];
  76. ]

  77. Print["x=" MatrixForm[s.a]];
  78. Print["Maxr="];
  79. maxr
  80. Print["Det[kmc]="];
  81. Det[kmc]
  82. Print["da"];
  83. da
复制代码

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2010-5-27 20:47 | 显示全部楼层

回复 楼主 如枫 的帖子

楼主能不能给一个具体一点的增量谐波平衡法的定义呀?
 楼主| 发表于 2010-6-2 08:32 | 显示全部楼层
楼上可以看看陈树辉的《强非线性振动系统的定量分析方法》,编程倒问题不大,关键是选初始解太难了,
我请教过陈树辉的学生,也说没什么好的办法!

评分

1

查看全部评分

发表于 2010-6-3 21:23 | 显示全部楼层

回复 板凳 如枫 的帖子

兄弟,你说的陈树辉《强非线性振动系统的定量分析方法》是电子版的吗?
 楼主| 发表于 2010-6-5 14:19 | 显示全部楼层
有纸质版的,没见到过电子版的
发表于 2010-8-28 18:39 | 显示全部楼层
收敛怎么解决?
发表于 2010-8-30 09:10 | 显示全部楼层
回复 如枫 的帖子

程序一般没什么问题,主要就是取值,最近我也在做,用的是多尺度法,也是总仿不出来,我觉得可能还是没取对值
   
发表于 2010-9-23 10:22 | 显示全部楼层
回复 kangarooli 的帖子

请问,你问题解决了吗?如何选取初值,有什么技巧没有呀?
发表于 2010-9-23 10:47 | 显示全部楼层
回复 liuwg14 的帖子

我的也没解决,也没什么技巧,我就是多取几个试的
发表于 2010-9-25 16:00 | 显示全部楼层
谢谢,我发现初值太敏感了.稍微有点改动,就得不到解. 我求解的是Mathieu方程!
发表于 2010-9-26 10:08 | 显示全部楼层
这个问题不解决,后期工作不好弄呀
 楼主| 发表于 2010-9-26 20:52 | 显示全部楼层
大家讨论一下初值的问题啊!能不能改进一下饭算法,使程序对初值不要这么敏感。
感觉数值解对初值选取有点帮助!
发表于 2010-9-26 21:53 | 显示全部楼层

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

个人认为还是从寻找初值的范围入手更加合理点
 楼主| 发表于 2010-10-4 20:15 | 显示全部楼层
自治系统好办一点,非自治系统的程序难写,总是不收敛!有好多文献用IHB法求非自治系统的解,不知道具体怎么操作的!
发表于 2010-10-9 16:32 | 显示全部楼层
如枫 发表于 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. epsilon=0.125;
  6. alpha=1.0;
  7. omega=2.0;
  8. gamma=1.0;
  9. beta=4.6;

  10. x=a0+a1*cos(t)+a2*cos(2*t)+b1*sin(t)+b2*sin(2*t);
  11. %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);

  12. % f=diff(diff(x,t),t)+2*epsilon*diff(x,t)-(alpha+beta*sin(omega*t))*x-gamma*x^3;
  13. %f=omega^2*diff(diff(x,t),t)+2*omega*epsilon*diff(x,t)-(alpha+beta*sin(t))*x-gamma*x^3
  14. d_x=diff(x,'t');
  15. dd_x=diff(d_x,'t');

  16. f=omega^2*dd_x+2*epsilon*omega*d_x-(alpha+beta*sin(t))*x+gamma*x^3

  17. f_1=simplify(int(f,t,0,2*pi)/(2*pi));
  18. f_cost=simplify(int(f*cos(t),t,0,2*pi)/pi);
  19. f_cos2t=simplify(int(f*cos(2*t),t,0,2*pi)/pi);
  20. f_sint=simplify(int(f*sin(t),t,0,2*pi)/pi);
  21. f_sin2t=simplify(int(f*sin(2*t),t,0,2*pi)/pi);

  22. % fun=[f_1;f_cost;f_cos2t;f_sint;f_sin2t]
  23. % x0(5)=0
  24. % Newtons(fun,x0)
复制代码
回复 支持 1 反对 0

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-3-28 18:52 , Processed in 0.166030 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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