简介：基于 Newmark β 方法的 MATLAB 实现，用于迭代求解双线性单自由度（SDOF）结构动力响应，适合本科、硕士阶段开展结构动力学与非线性地震响应分析教学或科研。压缩包共 3 个文件，含核心 .m 源码、运行结果对比图 .png 及输出数据 .csv，可直接在 MATLAB R2019a 中运行复现，便于对照检查每一步迭代结果。已有 121 人学习下载，代码体量小（18KB），结构清晰，特别适合初学者理解双线性滞回模型与数值积分方法的内在逻辑。借助这份资源，可以掌握 Newmark β 法的参数设置、双线性刚度切换处理以及迭代收敛判定技巧，同时利用 CSV 结果自行绘制响应曲线或扩展为多自由度模型，是一份简洁实用的教研参考。

1. 双线性 SDOF 的 Newmark-β 迭代求解：一份能跑通的地震时程分析代码

做结构抗震分析时，第一个非线性模型往往是双线性单自由度体系。用 Newmark-β 法做时程积分本来是标准操作，可一旦刚度进入屈服段，每一步都要迭代求解平衡方程，很多人的 MATLAB 代码就在这个环节开始乱跳、不收敛、滞回曲线变成一堆麻点。这份资源把整个求解器封装成可直接运行的 MATLAB 函数，适合刚接触非线性时程分析的研究生，也适合需要在小型工程里快速验证的工程师。它能解决的是：给你一套不需要从头攒的、带迭代收敛控制的 Newmark-β 双线性 SDOF 求解代码，同时让你能清楚看到每个参数对结果的影响。如果你正在被“位移漂移”或“迭代炸掉”折磨，这份代码值得花一晚上拆一遍。

2. Newmark-β 积分器与双线性本构：先把参数含义吃透再动手

拿到代码先别急着跑，Newmark-β 不是简单地套公式循环，两个积分参数和恢复力模型的配合直接决定你迭代是否收敛。这一章把运动方程、积分参数、双线性刚度三件事串起来，后面改参数时才不会瞎试。

2.1 积分参数 β、γ 的选型：平均加速度法为什么是默认选项

Newmark 法的核心是在一个时间步 Δt 内对加速度的分布做假设，然后推导出位移和速度的更新公式。最常见的写法是：

u_{n+1} = u_n + Δt · v_n + Δt² · [(1/2 − β) · a_n + β · a_{n+1}]

v_{n+1} = v_n + Δt · [(1 − γ) · a_n + γ · a_{n+1}]

β 和 γ 控制了数值积分的行为。γ 通常取 1/2，这时候积分不会引入人工阻尼；β 决定加速度在步内的插值方式。β=1/4 对应平均加速度法，积分无条件稳定，也就是步长取多大都不会出现数值发散，所以默认用它。β=1/6 对应线性加速度法，只有满足 Δt ≤ 0.551T 时才稳定，遇到周期很短的结构就容易翻车。

在双线性 SDOF 里，屈服后刚度变小，结构周期变长，平均加速度法的稳定性优势就更重要。我一般会把 β=1/4、γ=1/2 写死，只在对比精度时才去改 β。另一个容易被忽略的点是：γ 取 1/2 时没有数值阻尼，所以高频噪声不会被自动滤掉，这是后面避坑章要解决的问题。

提示：如果你发现位移时程尾部出现高频锯齿，先怀疑 β 和 γ 的取值，再怀疑刚度切换逻辑。

2.2 双线性恢复力模型：刚度突变点与屈服后刚度比

双线性模型的参数不多：初始弹性刚度 k₁、屈服位移 u_y、屈服后刚度 k₂ = α · k₁，其中 α 是屈服后刚度折减系数，常见取 0.02～0.1。恢复力 f_s 与位移 u 的关系不是单值函数，它取决于加载历史：从零开始加载，位移小于 u_y 时刚度是 k₁；超过 u_y 后刚度切到 k₂；卸载时刚度通常回到 k₁，形成滞回环。

这个模型难在“状态切换”。同一位移下，加载和卸载对应的恢复力不同，必须用一个状态变量记录当前是处于弹性加载、屈服加载、弹性卸载还是反向屈服。很多迭代发散就发生在刚度切换的那一步：你用了上一时刻的刚度来预测本步位移，结果预测值越过了屈服点，实际刚度和预测刚度不一致，Newton-Raphson 循环就得多跑几次。

在双线性 SDOF 中，切线刚度只可能是 k₁ 或 k₂，判断条件也简单：如果当前处于加载路径且位移超过历史最大位移，就用 k₂；否则用 k₁。不过卸载后再反向加载，如果反向位移也超过了反向屈服点，刚度同样要切成 k₂。所以代码里至少需要记录正向屈服位移、反向屈服位移、当前恢复力方向这三个量。

2.3 运动方程离散：从增量平衡方程到 Newton-Raphson 迭代

单自由度体系在任意时刻 t_{n+1} 必须满足动力平衡方程：

m · a_{n+1} + c · v_{n+1} + f_s(u_{n+1}) = F_{n+1}

其中 f_s 是恢复力，由双线性本构算出，不再是简单的 k·u。因为 f_s 是非线性的，所以这个方程不能直接解出 u_{n+1}，必须迭代。把 Newmark 的位移和速度更新公式代入，可以把 a_{n+1} 和 v_{n+1} 都写成 u_{n+1} 的函数，得到关于 u_{n+1} 的等效静力方程。

对每一步迭代，需要计算切线刚度 K_T（当前刚度的斜率，即 k₁ 或 k₂），再组装等效刚度：

K_eff = K_T + (γ/(β·Δt))·c + (1/(β·Δt²))·m

然后计算不平衡力 R = F_{n+1} − m·a_{n+1} − c·v_{n+1} − f_s(u_{n+1})，解出位移增量 Δu = R / K_eff，更新 u 后收敛判断。这套流程就是 Newton-Raphson 法在时程分析里的标准做法。

理解到这里，你再看压缩包里的代码就不会晕：主程序负责时间循环，子函数负责每一步的迭代和刚度更新。接下来一章直接拆解这份 MATLAB 代码的骨架。

3. MATLAB 实现：从时程循环到滞回曲线的可复现代码

这一章是全文动手的核心。我不会把完整上千行贴出来，而是抽出三个关键代码块：主循环骨架、迭代求解函数、后处理片段。你拿到压缩包后可以用这些代码块对照原文件，看明白每段在干什么，再改参数。

3.1 主程序骨架：状态变量初始化与时间循环

% newmark_bilinear_sdof.m % 双线性SDOF Newmark-beta迭代求解，主程序 clc; clear; close all; % --- 结构参数 --- m = 1.0; % 质量，单位自定义 k1 = 100.0; % 初始弹性刚度 alpha = 0.05; % 屈服后刚度比 uy = 0.02; % 屈服位移 zeta = 0.05; % 阻尼比 c = 2*zeta*sqrt(m*k1);% 粘性阻尼系数（基于初始刚度） % --- 地震波参数 --- dt = 0.01; % 时间步长，秒 ag = load('elcentro.txt'); % 地震加速度时程，单位与结构单位一致 N = length(ag); % --- 状态变量初始化 --- u = zeros(N,1); % 位移 v = zeros(N,1); % 速度 a = zeros(N,1); % 加速度 fs = zeros(N,1); % 恢复力 hist_state = zeros(N,1); % 记录刚度状态，1为弹性/卸载，2为屈服加载 % --- 时间循环 --- for i = 2:N % 当前时刻激励 F = -m * ag(i); % 地震作用下等效惯性力 % 调用单步迭代求解器 [u(i), v(i), a(i), fs(i), hist_state(i)] = ... newmark_step(m, c, k1, alpha, uy, u(i-1), v(i-1), a(i-1), ... F, dt, hist_state(i-1)); end

逻辑说明：主程序把结构参数拆开放在开头，避免在函数内部反复传参。注意阻尼系数 c 用初始刚度计算，这在双线性问题里是常规做法。地震波加载ag时，必须保证物理单位一致，例如质量取 kg、位移取 m 时，加速度单位应是 m/s²。循环体内每一结构步调用newmark_step，把上一时刻的位移、速度、加速度和刚度状态传进去，返回本时刻结果。

参数说明：alpha越小屈服后刚度越软，滞回环越“胖”，但迭代越容易震荡；uy越小结构越早屈服，非线性程度越高。dt=0.01只适合周期在 0.2s 以上的结构，后面第 4 章会专门讲步长上限。

3.2 迭代求解器函数：刚度更新与不平衡力修正

function [u_new, v_new, a_new, fs_new, state_new] = ... newmark_step(m, c, k1, alpha, uy, u_old, v_old, a_old, F, dt, state_old) % Newmark-beta单步迭代，采用Newton-Raphson修正 beta = 1/4; gamma = 1/2; % 常数 a0 = 1/(beta*dt^2); a1 = gamma/(beta*dt); a2 = 1/(beta*dt); a3 = 1/(2*beta) - 1; % 预测位移（先用上一步位移作为迭代初值） u_iter = u_old; % 恢复力与切线刚度初始计算 kt = k1; fs_iter = bilinear_force(u_iter, k1, alpha, uy, state_old); if abs(u_iter) > uy kt = alpha * k1; % 简化判断，实际需结合状态 end % Newton-Raphson迭代 for j = 1:20 % 根据当前u_iter计算a和v a_iter = a0*(u_iter - u_old) - a2*v_old - a3*a_old; v_iter = v_old + (1-gamma)*dt*a_old + gamma*dt*a_iter; % 不平衡力 R = F - m*a_iter - c*v_iter - fs_iter; % 等效刚度 K_eff = kt + a1*c + a0*m; % 位移修正量 du = R / K_eff; u_iter = u_iter + du; % 更新恢复力和切线刚度 [fs_iter, kt, state_new] = bilinear_force(u_iter, k1, alpha, uy, state_old); % 收敛判断（同时检查位移修正量和残余力） if abs(du) < 1e-8 * max(1, abs(u_iter)) && abs(R) < 1e-8 * max(1, abs(F)) break; end end % 输出本步结果 u_new = u_iter; a_new = a0*(u_new - u_old) - a2*v_old - a3*a_old; v_new = v_old + (1-gamma)*dt*a_old + gamma*dt*a_new; fs_new = fs_iter; end

逻辑说明：这个函数是单步迭代的核心。先用上一步位移做初值，然后进入最多 20 次 Newton-Raphson 循环。每次循环先根据当前位移通过 Newmark 关系式计算加速度和速度，再算不平衡力 R，解出修正位移量 du。重点在于收敛判据同时检查位移修正量和残余力，避免单看位移误差时残余力还很大。

参数说明：迭代上限 20 次看起来不大，但对双线性 SDOF 这种分段线性问题通常 3～5 次就收敛。如果达到 20 次还没收敛，程序没有报错而是直接输出结果，这样可能掩盖问题，你可以在循环外加一个if j == 20 && abs(R) > 1e-6 * max(1, abs(F))的警告，方便定位发散点。

3.3 后处理：提取滞回曲线、位移时程与能量指标

% 后处理脚本 figure(1); plot(t, u, 'b-', 'LineWidth', 1.2); xlabel('Time (s)'); ylabel('Displacement'); title('Displacement Time History'); grid on; figure(2); plot(u, fs, 'r-'); xlabel('Displacement'); ylabel('Restoring Force'); title('Hysteresis Loop'); grid on; % 计算累积滞回耗能 energy = zeros(N,1); for i = 2:N avg_fs = 0.5 * (fs(i) + fs(i-1)); du = u(i) - u(i-1); energy(i) = energy(i-1) + avg_fs * du; end

逻辑说明：第一条绘图命令画位移时程，重点看是否有不可恢复的漂移。第二条画滞回曲线，理想双线性滞回环应该是一个倾斜的平行四边形，如果出现蝴蝶状或锯齿状就要回头查迭代。能量累计用梯形公式，把每两步之间的恢复力和位移增量乘起来累加，这个量用于第 6 章的验证。

参数说明：这里energy是每条边积分，如果加载过程中有卸载再加载，滞回环面积就是进出的差值，能量曲线应单调增长。后处理不需要太花哨，能把这几个量画出来就足够定位 90% 的问题。

4. 参数配置与调试：收敛容差、阻尼和地震波怎么配合

代码能跑只是第一步，参数配不对，结果照样不能用。这一章给出我在实际使用中最常调整的几组参数，以及它们背后的理由。

4.1 收敛容差与最大迭代次数：工程上常用的经验范围

收敛容差不是越小越好。太小的容差会导致迭代在刚度突变点附近反复震荡走不出去，太大的容差又会把误差累积进位移时程。常见做法是位移修正量相对容差取 1e-8，残余力相对值取 1e-8。如果你发现某几步迭代次数突然超过 10 次，优先检查该时刻是不是刚发生刚度切换。

最大迭代次数在双线性 SDOF 里通常取 20 足够，因为解析求导的切线刚度非常精确。如果你把代码改成了其他复杂本构，比如考虑退化或捏缩，最大迭代次数建议加到 50，并且加入发散保护：if isnan(du) || isinf(du), error('迭代发散')。这不影响正常工况，但能在参数极端时救你一命。

4.2 质量、刚度与阻尼矩阵的组装：避免硬编码

这份资源针对 SDOF，只有一个质量 m、一个阻尼 c 和一个刚度 k，看起来简单，但单位问题常常让人翻车。很多新手会把质量直接写成重量，导致周期偏大。正确的做法是：加速度用 m/s²，力用 N，质量用 kg；如果用 cm 做位移单位，弹性模量和力的单位要配套换算。

阻尼是容易被忽略的变量。双线性问题中阻尼一般基于初始刚度 k₁ 计算，也就是 c = 2ζ√(m·k₁)，而不是用当前切线刚度。如果阻尼随刚度切换，每一步迭代的阻尼力会跳变，反而引入虚假能量。压缩包里的代码应该也是这样写的，你可以验证一下。

注意：如果你输入的地震波本身已经做了单位换算，千万不能再乘重力加速度 g。我见过有人把 El Centro 波的单位误当成 m/s²，直接放进去，结果位移被放大了近十倍。

4.3 时程步长的选择：与地震波采样率和周期匹配

时间步长 dt 直接决定数值色散。对于 Newmark 平均加速度法，建议 dt ≤ T_min / 10，其中 T_min 是结构在整个非线性过程中可能出现的最短周期。双线性模型屈服后刚度变小，周期变大，所以真正限制步长的是弹性段周期。如果结构周期是 0.2s，dt 取 0.01s 是可以的；如果周期只有 0.05s，dt 至少应加密到 0.002s。

地震波采样率也要检查。常见的地震波文件采样间隔是 0.005s 或 0.01s，如果你的积分步长和采样间隔不一致，需要在加载时做线性插值。代码里如果直接load然后用原始点列当时间序列，大概率会忽略插值，这是导致激励与结构响应错位的隐藏坑。下表是常见参数范围：

5. 避坑指南：Newmark 迭代求解里最常见的五个翻车现场

这一章是我自己在调代码时踩过最多次的坑，每一条都是“现象 → 原因 → 解决”的结构。你遇到类似症状时可以直接对照排查。

5.1 残余力不归零导致位移漂移

现象：位移时程在强震结束后不回零，而是随着时间越来越大，形成明显的基线漂移。

原因：Newton-Raphson 迭代跳出循环时，残余力还比较大，导致当前步的加速度偏大，速度增量也偏大。这些误差在后续每个时间步里不断累积，尤其当结构进入塑性后，残余力方向与位移方向一致时，漂移会被放大。

解决：在收敛判断里同时加入残余力绝对值和相对值判断。相对值建议用abs(R) < 1e-8 * max(1, abs(F))而不是只比位移增量。另外，把迭代收敛后的恢复力fs_new重新赋值为平衡力F - m*a - c*v，这样能强制恢复力满足离散方程，避免滞回曲线不闭合。

5.2 刚度切换瞬间迭代发散

现象：迭代次数在某一步突然达到上限，或者位移增量出现异常大值。

原因：双线性刚度在屈服点不可微，用初值 u_old 计算切线刚度时，如果 u_old 在弹性侧但预测位移越过屈服点，Newton-Raphson 迭代会先用弹性刚度走一大步，再切回屈服刚度，来回震荡。这种分段线性问题在刚度突变点附近容易进入死循环。

解决：在迭代循环内每计算一次恢复力都要重新判断当前切线刚度，而不是在循环外算一次就固定住。另外，可以在迭代每步限制位移增量幅值，比如du = max(min(du, 10*uy), -10*uy)，防止单步跳出本构域。这种做法叫“增量限定”，对保持迭代稳定很有效。

5.3 β 取值不当导致高频响应失真

现象：位移时程与解析解对比时，峰值一致但曲线明显更毛糙，甚至出现虚假的高频波动。

原因：有人为了“更精确”把 β 改成 0，变成显式中心差分法，或者把 γ 从 1/2 改成其他值引入人工阻尼。Newmark 法的稳定性与 β 直接相关，β 过小时算法变成有条件稳定，步长一大就跑出高频假响应。

解决：双线性 SDOF 老老实实用 β=1/4，γ=1/2。如果你需要高频滤波，正确的做法是在地震波激励上做低通滤波，而不是靠积分器引入人工阻尼。人工阻尼会削弱真实地震响应，尤其在峰值加速度附近。

5.4 收敛判据选位移还是力

现象：程序汇报“已收敛”，但滞回曲线在屈服点附近呈现一个尖角而不是平滑过渡。

原因：只检查位移修正量小于容差就跳出循环，但此时残余力仍然很大。位移增量小并不代表平衡方程满足，因为刚度为双线性时，很小的位移误差也会引起较大的力误差。

解决：位移和残余力双准则，哪个不满足都要继续迭代。我一般写if abs(du) < 1e-8 && abs(R) < 1e-8 * max(1, abs(F))，两个条件都要满足才退出。这样会增加几次迭代，但对结果质量的提升非常明显。

5.5 与 OpenSees 或其他软件结果对不上时先查单位

现象：同样的结构参数和地震波，MATLAB 算出的最大位移是 OpenSees 的十分之一或十倍，滞回环形状类似但幅值差很多。

原因：单位不一致。OpenSees 用 N 和 m 时，弹性模量单位是 Pa；而很多同学在 MATLAB 里用 kN 和 mm 混算，导致质量、刚度、加速度三个量互相矛盾。更隐蔽的是地震波单位，有的文件写的是 “gal”（cm/s²），如果不除以 100 换算成 m/s²，结果就会大两三个量级。

解决：核对五件事：质量单位、位移单位、时间单位、力单位、加速度单位。建议把所有物理量都标准化为 m、kg、s、N，地震波加速度统一换算成 m/s²。压缩包代码里如果地震波读取没有单位换算，你跑之前先检查文件头的加速度范围，El Centro 的典型峰值约 0.3g 也就是约 3 m/s²，如果你的波峰峰值是 300 左右，那单位基本是 gal。

6. 进阶验证：用刚度折减后的周期和滞回环面积给代码把关

代码调通之后，建议做两个验证性实验，确认你的求解器没在“看起来正常”的表面下藏误差。我在拿到这份资源时做的第一件事就是把屈服位移设成极大，让结构全程弹性，然后用解析解对比，这样能先把积分器本身的问题排除掉。

% 验证1：弹性验证 uy = 1e10; % 永远不屈服 alpha = 0; % 屈服后刚度无意义 % 运行求解器后，提取最大位移 % 解析解：单自由度正弦激励下的稳态响应

令屈服位移极大时，双线性模型退化为弹性，积分器的精度可以直接用谐波激励的稳态响应检验。如果最大位移误差超过 1%，先怀疑阻尼比 ζ 或 β 取错，再去查恢复力函数是否意外触发了屈服逻辑。

第二个验证是滞回环面积。双线性 SDOF 在一个完整滞回周期内，理论耗能等于滞回环包围的面积，也就是 4 倍的“屈服位移×屈服后力差”再修正。你可以手动加一个正弦循环加载，把求解器输出的 u 和 fs 画出来后用polyarea(u, fs)计算面积，再和手算值对比。误差应小于迭代容差的影响范围。这个验证能直接暴露状态切换逻辑的问题，如果面积偏大或偏小，说明卸载刚度或者反向屈服判断有误。

我自己的习惯是每次跑完一个工况，强制看三个量：第一个是最大位移，第二个是滞回曲线是否闭合，第三个是能量曲线是否单调递增。这三个量任何一个异常，我都不会再往下做后处理，而是回头查参数。从那以后，我每次拿到新的 SDOF 求解代码，都会先跑两个半小时的验证实验，而不是直接丢进论文里算结果。这份压缩包里的代码也是一样，你下载后花一晚上跑通和验证，后面做参数分析和参数识别时就再也不用担心求解器本身是黑匣子了。希望帮到你。

本文还有配套的精品资源，点击获取