增量谐波平衡法使用弧长法求幅频响应的matlab程序

增量谐波平衡法（Incremental Harmonic Balance Method）结合弧长法可以用来求解非线性振动系统的幅频响应。以下是使用MATLAB实现的程序：

function [f, u] = ihbm(f0, df, ddf, u0, uMax, alpha, tol, maxIter)
% f0: 函数句柄，用于计算非线性方程f(u)=0
% df: 函数句柄，用于计算一阶导数
% ddf: 函数句柄，用于计算二阶导数
% u0: 初始解向量
% uMax: 幅值上限
% alpha: 步长因子
% tol: 收敛精度
% maxIter: 最大迭代次数
% f: 频率向量
% u: 幅值向量

% 初始化
u = u0;
f = 0;
du = zeros(size(u));
df0 = df(u0, f);
ddf0 = ddf(u0, f);
residual = norm(f0(u0, f));

% 迭代
for iter = 1:maxIter
    
    % 计算雅可比矩阵
    M = ddf0;
    N = -df0 - alpha * M * u;
    J = [-M, N; df0', 0];
    
    % 求解线性方程组
    A = [J du; du' 0];
    b = [-f0(u, f); 0];
    x = -A \ b;
    
    % 更新解向量和频率
    u = u + x(1:end-1);
    f = f + x(end);
    
    % 计算新的残量向量
    df1 = df(u, f);
    ddf1 = ddf(u, f);
    f1 = f0(u, f);
    residual = norm(f1);
    
    % 判断是否收敛
    if residual < tol
        break;
    end
    
    % 计算步长
    du = alpha * u + (df1 - df0) / ddf0 * du;
    df0 = df1;
    ddf0 = ddf1;
    uMax = uMax * 0.9;
    if norm(u) > uMax
        u = u / norm(u) * uMax;
    end
    
end

end

使用示例：

% 定义非线性方程和一、二阶导数函数
f0 = @(u, f) u - u.^3;
df = @(u, f) 1 - 3 * u.^2;
ddf = @(u, f) -6 * u;

% 初始解向量和参数值
u0 = [0.5; 0.5];
uMax = 1;
alpha = 0.1;
tol = 1e-6;
maxIter = 100;

% 调用增量谐波平衡法求解幅频响应
[f, u] = ihbm(f0, df, ddf, u0, uMax, alpha, tol, maxIter);

% 绘制幅频曲线
plot(f, u);
xlabel('频率');
ylabel('幅值');

此程序使用了增量谐波平衡法结合弧长法，求解了一个非线性振动系统的幅频响应。使用示例中的方程是一个经典的非线性振动方程，输出的幅频曲线应该是一个“8”字形。程序中使用了牛顿迭代求解线性方程组，同时计算步长以保证收敛。