关于增量谐波平衡法（IHB）的编程问题
有没有人编过增量谐波平衡法的Mathematica程序？我最近在编程序，一直都不行，想请教一下下面是我编的Van der pol方程程序，但是好像不收敛，不知道是程序有问题还是初值选取不好Unprotect["`*"];
ClearAll["`*"];
(*输入初始数据*)
nc = 1; ns = 1; n = nc + ns;
omega0 = 0.94;
epsilon = 1.0;
(*输入初始矩阵*)
cba = epsilon; (*输入C拔矩阵*)
kba = 1; (*输入K拔矩阵*)
a = Array; (*输入A的初始矩阵*)
a[] = 2;
Do] = 0.241, {i, 1, nc}];
Do] = 0.126, {i, nc + 2, n}];
(*主程序*)
s = Array;
Do = Cos, {i, nc}];
Do = Sin[(i - nc)*t], {i, nc + 1, n}];
s1 = Dt;
s2 = Dt;
st = Transpose;
m0 = st.s2;
m = Integrate;
c0 = cba st.s1;
c = Integrate;
k0 = kba st.s;
k = Integrate;
r = Array;
Do = 1, {i, n}];
maxr = 1;
nn = 1;
x = s.a;
c3ba = epsilon x^2;
x1 = Dt;
k2ba = epsilon x x1;
c30 = c3ba[] st.s1;
k20 = k2ba[] st.s;
c3 = Integrate;
k2 = Integrate;
kmc = omega0^2 m + omega0 (c3 - c) + k + 2 omega0 k2;
rmc = -(2 omega0 m + c3 - c).a;
r = -(omega0^2 m + omega0 (c3 - c) + k + k2).a;
While[maxr > 0.00001,
(*将kmc中第一列元素换为rmc*)
Do] = -rmc[], {i, 1, n}];
(*定义方程*)
(*left=kmc.da;
right=r-domega*rmc;*)
(*If<0.1*10^(-10),
omega0=omega0+domega(*如果omega0增大不行，改成减小*);Continue[]];*)
da = LinearSolve] &)];
(*Print["a0="]
a0
a0+da[]
omega0
omega0+domega*)
omega0 = omega0 + da[];
Do] = a[] + da[], {i, 1, nc}];
Do] = a[] + da[], {i, nc + 2, n}];
(*omega0=omega0+domega;*)
(*a0
omega0*)
x = s.a;
c3ba = epsilon x^2;
x1 = Dt;
k2ba = epsilon x x1;
c30 = c3ba[] st.s1;
k20 = k2ba[] st.s;
c3 = Integrate;
k2 = Integrate;
kmc = omega0^2 m + omega0 (c3 - c) + k + 2 omega0 k2;
rmc = -(2 omega0 m + c3 - c).a;
r = -(omega0^2 m + omega0 (c3 - c) + k + k2).a;
maxr = Max]];
nn = nn + 1;
If; Break[]];
]
Print["x=" MatrixForm];
Print["Maxr="];
maxr
Print["Det="];
Det
Print["da"];
da
回复 楼主 如枫 的帖子
楼主能不能给一个具体一点的增量谐波平衡法的定义呀？ 楼上可以看看陈树辉的《强非线性振动系统的定量分析方法》，编程倒问题不大，关键是选初始解太难了，我请教过陈树辉的学生，也说没什么好的办法！
回复 板凳 如枫 的帖子
兄弟，你说的陈树辉《强非线性振动系统的定量分析方法》是电子版的吗？ 有纸质版的，没见到过电子版的 收敛怎么解决？ 回复 如枫 的帖子程序一般没什么问题，主要就是取值，最近我也在做，用的是多尺度法，也是总仿不出来，我觉得可能还是没取对值
回复 kangarooli 的帖子
请问,你问题解决了吗?如何选取初值,有什么技巧没有呀? 回复 liuwg14 的帖子
我的也没解决，也没什么技巧，我就是多取几个试的 谢谢,我发现初值太敏感了.稍微有点改动,就得不到解. 我求解的是Mathieu方程! {:{44}:}这个问题不解决,后期工作不好弄呀 大家讨论一下初值的问题啊！能不能改进一下饭算法，使程序对初值不要这么敏感。
感觉数值解对初值选取有点帮助！ 如枫 发表于 2010-9-26 20:52 static/image/common/back.gif
大家讨论一下初值的问题啊！能不能改进一下饭算法，使程序对初值不要这么敏感。
感觉数值解对初值选取有点 ...
这个问题估计不好办，对初值是否敏感是有方程自身的特性决定的
如果改变了那还是原来的方程吗？
个人认为还是从寻找初值的范围入手更加合理点 自治系统好办一点，非自治系统的程序难写，总是不收敛！有好多文献用IHB法求非自治系统的解，不知道具体怎么操作的！ 如枫 发表于 2010-10-4 20:15 static/image/common/back.gif
自治系统好办一点，非自治系统的程序难写，总是不收敛！有好多文献用IHB法求非自治系统的解，不知道具体怎么 ...
我这儿有一个matlab版的IHB法程序，网站发现的，下载了还没用过，贴出来大家看看是否好使function y=func()
clc
clear
syms t a0 a1 a2 a3 a4 a5 a6 a7 a8 b1 b2 b3 b4 b5 b6 b7 b8;
epsilon=0.125;
alpha=1.0;
omega=2.0;
gamma=1.0;
beta=4.6;
x=a0+a1*cos(t)+a2*cos(2*t)+b1*sin(t)+b2*sin(2*t);
%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);
% f=diff(diff(x,t),t)+2*epsilon*diff(x,t)-(alpha+beta*sin(omega*t))*x-gamma*x^3;
%f=omega^2*diff(diff(x,t),t)+2*omega*epsilon*diff(x,t)-(alpha+beta*sin(t))*x-gamma*x^3
d_x=diff(x,'t');
dd_x=diff(d_x,'t');
f=omega^2*dd_x+2*epsilon*omega*d_x-(alpha+beta*sin(t))*x+gamma*x^3
f_1=simplify(int(f,t,0,2*pi)/(2*pi));
f_cost=simplify(int(f*cos(t),t,0,2*pi)/pi);
f_cos2t=simplify(int(f*cos(2*t),t,0,2*pi)/pi);
f_sint=simplify(int(f*sin(t),t,0,2*pi)/pi);
f_sin2t=simplify(int(f*sin(2*t),t,0,2*pi)/pi);
% fun=
% x0(5)=0
% Newtons(fun,x0)
页:
[1]
2