MATLAB 常微分方程数值求解算法探索：以两自由度无阻尼振动系统为例

MATLAB常微分方程数值求解算法程序龙格库塔法、威尔逊法、纽马克法、中心差分法以两自由度无阻尼振动系统为例在MATLAB中建模并编制数值计算输出四种算法下物块的位移、速度和加速度曲线后续可在此基础上继续开展算法求解误差对比。在动力学分析领域常微分方程ODE的数值求解是一项关键任务。今天咱们就来聊聊 MATLAB 中针对常微分方程的几种常见数值求解算法包括龙格 - 库塔法、威尔逊法、纽马克法以及中心差分法。咱们以两自由度无阻尼振动系统为实例看看如何在 MATLAB 里建模并用这几种算法算出物块的位移、速度和加速度曲线后续还能对这些算法的求解误差做个对比。两自由度无阻尼振动系统的数学模型两自由度无阻尼振动系统可以用下面的二阶常微分方程组来描述\[m1\ddot{x}1 k1 x1 - k2 (x2 - x_1) 0\]\[m2\ddot{x}2 k2 (x2 - x_1) 0\]为了能在 MATLAB 里求解我们把它转化为一阶常微分方程组。设 \( y1 x1 \)\( y2 \dot{x}1 \)\( y3 x2 \)\( y4 \dot{x}2 \)则方程组变为\[\dot{y}1 y2MATLAB常微分方程数值求解算法程序龙格库塔法、威尔逊法、纽马克法、中心差分法以两自由度无阻尼振动系统为例在MATLAB中建模并编制数值计算输出四种算法下物块的位移、速度和加速度曲线后续可在此基础上继续开展算法求解误差对比。\]\[\dot{y}2 \frac{- (k1 k2)y1 k2 y3}{m_1}\]\[\dot{y}3 y4\]\[\dot{y}4 \frac{k2 y1 - k2 y3}{m2}\]龙格 - 库塔法龙格 - 库塔法是一种常用且高精度的数值求解方法。在 MATLAB 里我们可以利用ode45函数来实现它就是基于龙格 - 库塔法的。% 参数设置 m1 1; m2 1; k1 100; k2 100; tspan 0:0.01:10; % 时间范围 y0 [0; 0; 0; 0]; % 初始条件 % 定义微分方程 odefun (t,y) [y(2); (- (k1 k2)*y(1) k2*y(3))/m1; y(4); (k2*y(1) - k2*y(3))/m2]; % 使用ode45求解 [t,y] ode45(odefun,tspan,y0); % 提取位移、速度和加速度 x1 y(:,1); v1 y(:,2); a1 (- (k1 k2)*y(:,1) k2*y(:,3))./m1; x2 y(:,3); v2 y(:,4); a2 (k2*y(:,1) - k2*y(:,3))./m2; % 绘制曲线 figure; subplot(3,1,1); plot(t,x1,t,x2); title(位移曲线); xlabel(时间 t (s)); ylabel(位移 x (m)); legend(物块1位移,物块2位移); subplot(3,1,2); plot(t,v1,t,v2); title(速度曲线); xlabel(时间 t (s)); ylabel(速度 v (m/s)); legend(物块1速度,物块2速度); subplot(3,1,3); plot(t,a1,t,a2); title(加速度曲线); xlabel(时间 t (s)); ylabel(加速度 a (m/s^2)); legend(物块1加速度,物块2加速度);代码分析首先设定了系统参数 \( m1 \)、\( m2 \)、\( k1 \)、\( k2 \)时间范围tspan和初始条件y0。然后定义了微分方程odefun这就是前面转化后的一阶常微分方程组。接着使用ode45函数求解微分方程得到时间t和状态变量y。从y中提取出物块 1 和物块 2 的位移、速度和加速度。最后用subplot函数绘制出位移、速度和加速度曲线。威尔逊法威尔逊法是一种逐步积分法常用于结构动力学分析。以下是威尔逊法求解两自由度无阻尼振动系统的代码示例% 参数设置 m1 1; m2 1; k1 100; k2 100; tspan 0:0.01:10; % 时间范围 dt tspan(2)-tspan(1); % 时间步长 y0 [0; 0; 0; 0]; % 初始条件 % 初始化变量 n length(tspan); x1 zeros(n,1); v1 zeros(n,1); a1 zeros(n,1); x2 zeros(n,1); v2 zeros(n,1); a2 zeros(n,1); x1(1) y0(1); v1(1) y0(2); a1(1) (- (k1 k2)*y0(1) k2*y0(3))/m1; x2(1) y0(3); v2(1) y0(4); a2(1) (k2*y0(1) - k2*y0(3))/m2; % 威尔逊法迭代 for i 1:n-1 % 预测 a1_pred a1(i); a2_pred a2(i); v1_pred v1(i) dt*a1(i); v2_pred v2(i) dt*a2(i); x1_pred x1(i) dt*v1(i) 0.5*dt^2*a1(i); x2_pred x2(i) dt*v2(i) 0.5*dt^2*a2(i); % 校正 M [m1 0 0 0; 0 m1 0 0; 0 0 m2 0; 0 0 0 m2]; K [k1 k2 -k2 0 0; -k2 k1 k2 0 0; 0 0 k2 -k2; 0 0 -k2 k2]; F [0; 0; 0; 0]; a M \ (F - K * [x1_pred; x2_pred]); a1(i1) a(1); a2(i1) a(3); v1(i1) v1_pred 0.5*dt*(a1(i1) - a1_pred); v2(i1) v2_pred 0.5*dt*(a2(i1) - a2_pred); x1(i1) x1_pred dt*v1_pred 0.1666667*dt^2*(a1(i1) 2*a1_pred); x2(i1) x2_pred dt*v2_pred 0.1666667*dt^2*(a2(i1) 2*a2_pred); end % 绘制曲线 figure; subplot(3,1,1); plot(tspan,x1,tspan,x2); title(位移曲线); xlabel(时间 t (s)); ylabel(位移 x (m)); legend(物块1位移,物块2位移); subplot(3,1,2); plot(tspan,v1,tspan,v2); title(速度曲线); xlabel(时间 t (s)); ylabel(速度 v (m/s)); legend(物块1速度,物块2速度); subplot(3,1,3); plot(tspan,a1,tspan,a2); title(加速度曲线); xlabel(时间 t (s)); ylabel(加速度 a (m/s^2)); legend(物块1加速度,物块2加速度);代码分析同样先设定参数、时间范围和初始条件。初始化位移、速度和加速度变量并根据初始条件算出初始值。在迭代循环中先进行预测步骤根据上一步的加速度、速度和位移预测下一步的值。然后通过系统的质量矩阵M、刚度矩阵K和外力向量F进行校正得到更准确的加速度、速度和位移。最后绘制曲线展示结果。纽马克法纽马克法也是一种逐步积分方法下面是它的实现代码% 参数设置 m1 1; m2 1; k1 100; k2 100; tspan 0:0.01:10; % 时间范围 dt tspan(2)-tspan(1); % 时间步长 y0 [0; 0; 0; 0]; % 初始条件 % 初始化变量 n length(tspan); x1 zeros(n,1); v1 zeros(n,1); a1 zeros(n,1); x2 zeros(n,1); v2 zeros(n,1); a2 zeros(n,1); x1(1) y0(1); v1(1) y0(2); a1(1) (- (k1 k2)*y0(1) k2*y0(3))/m1; x2(1) y0(3); v2(1) y0(4); a2(1) (k2*y0(1) - k2*y0(3))/m2; % 纽马克法参数 beta 0.25; gamma 0.5; % 纽马克法迭代 for i 1:n-1 % 计算中间变量 A M gamma*dt*C beta*dt^2*K; b F M*(v1(i) (1 - gamma)*dt*a1(i)) C*(x1(i) dt*v1(i) (0.5 - beta)*dt^2*a1(i)); % 求解位移 x A \ b; x1(i1) x(1); x2(i1) x(3); % 计算速度和加速度 a1(i1) (x1(i1) - x1(i) - dt*v1(i)) / (beta*dt^2); a2(i1) (x2(i1) - x2(i) - dt*v2(i)) / (beta*dt^2); v1(i1) v1(i) (1 - gamma)*dt*a1(i) gamma*dt*a1(i1); v2(i1) v2(i) (1 - gamma)*dt*a2(i) gamma*dt*a2(i1); end % 绘制曲线 figure; subplot(3,1,1); plot(tspan,x1,tspan,x2); title(位移曲线); xlabel(时间 t (s)); ylabel(位移 x (m)); legend(物块1位移,物块2位移); subplot(3,1,2); plot(tspan,v1,tspan,v2); title(速度曲线); xlabel(时间 t (s)); ylabel(速度 v (m/s)); legend(物块1速度,物块2速度); subplot(3,1,3); plot(tspan,a1,tspan,a2); title(加速度曲线); xlabel(时间 t (s)); ylabel(加速度 a (m/s^2)); legend(物块1加速度,物块2加速度);代码分析设定参数、时间范围和初始条件与前面类似。初始化变量并根据初始条件算出初始值。定义纽马克法的参数beta和gamma。在迭代循环中先计算中间矩阵A和向量b然后求解位移。根据位移计算速度和加速度最后绘制曲线。中心差分法中心差分法是一种简单直观的数值求解方法。代码如下% 参数设置 m1 1; m2 1; k1 100; k2 100; tspan 0:0.01:10; % 时间范围 dt tspan(2)-tspan(1); % 时间步长 y0 [0; 0; 0; 0]; % 初始条件 % 初始化变量 n length(tspan); x1 zeros(n,1); v1 zeros(n,1); a1 zeros(n,1); x2 zeros(n,1); v2 zeros(n,1); a2 zeros(n,1); x1(1) y0(1); v1(1) y0(2); a1(1) (- (k1 k2)*y0(1) k2*y0(3))/m1; x2(1) y0(3); v2(1) y0(4); a2(1) (k2*y0(1) - k2*y0(3))/m2; % 中心差分法迭代 for i 2:n-1 % 计算加速度 a1(i) (- (k1 k2)*x1(i) k2*x2(i))/m1; a2(i) (k2*x1(i) - k2*x2(i))/m2; % 计算位移 x1(i1) 2*x1(i) - x1(i-1) dt^2*a1(i); x2(i1) 2*x2(i) - x2(i-1) dt^2*a2(i); % 计算速度 v1(i) (x1(i1) - x1(i-1)) / (2*dt); v2(i) (x2(i1) - x2(i-1)) / (2*dt); end % 绘制曲线 figure; subplot(3,1,1); plot(tspan,x1,tspan,x2); title(位移曲线); xlabel(时间 t (s)); ylabel(位移 x (m)); legend(物块1位移,物块2位移); subplot(3,1,2); plot(tspan,v1,tspan,v2); title(速度曲线); xlabel(时间 t (s)); ylabel(速度 v (m/s)); legend(物块1速度,物块2速度); subplot(3,1,3); plot(tspan,a1,tspan,a2); title(加速度曲线); xlabel(时间 t (s)); ylabel(加速度 a (m/s^2)); legend(物块1加速度,物块2加速度);代码分析同样设定参数、时间范围和初始条件并初始化变量。在迭代循环中先根据当前位移计算加速度。然后利用中心差分公式计算下一个时刻的位移。最后根据位移计算速度完成后绘制曲线。后续的算法求解误差对比通过上面的代码我们已经得到了四种算法下物块的位移、速度和加速度曲线。接下来就可以对比它们的求解误差了。一种常见的方法是与精确解对比如果精确解已知的话或者将一种高精度算法的结果作为参考解对比其他算法与参考解的误差。例如我们可以把龙格 - 库塔法ode45的结果作为参考解计算其他三种算法在每个时间点的误差然后绘制误差曲线这样就能直观地看出哪种算法在这个系统中的精度更高。具体实现这里就不再展开啦感兴趣的小伙伴可以自己动手试试。希望这篇博文能让大家对 MATLAB 中这几种常微分方程数值求解算法在两自由度无阻尼振动系统中的应用有更清晰的认识。大家要是有啥问题或者想法欢迎在评论区留言交流。