matlab代码求助！ - matlab吧

level 2

我需要求一个比较复杂的梁的振动方程，下面是这个问题的描述：

我查了很多论文，很简略地写了基于行列式|A|=0，就可以求出未知量α，但是都没有写出具体的求法。我在matlab里也尝试用类似的方法，把分块矩阵A写出来，然后令它的行列式等于零det（A）=0，然后通过fsolve求解，但是结果总是报错，请问各位大佬有没有知道类似问题该怎么在matlab里求解吗？

2025年01月02日 06点01分 1

吧务

level 12

编程的德彪

报错了啥

2025年01月02日 06点01分 2

仰望着丿星空

也不是报错吧，就是提示方程在初始点处解出。不知道是不是我的方程包含了太多的超越函数。

2025年01月02日 07点01分

编程的德彪

@仰望着丿星空是那种黄色的警告吗

2025年01月02日 08点01分

仰望着丿星空

@编程的德彪不是黄色的警告，我写的代码放在后面了，大佬方便帮我看一下吗，谢谢 [太开心]

2025年01月02日 08点01分

level 2

仰望着丿星空楼主

也不是报错吧，就是提示方程在初始点处解出。但是手算的话得出的并不是这个结果。
【Equation solved at initial point.
fsolve completed because the vector of function values at the initial point
is near zero as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>】

2025年01月02日 07点01分 3

level 2

仰望着丿星空楼主

顶

2025年01月02日 08点01分 4

level 2

仰望着丿星空楼主

clear all; clc;
syms omega;
% 梁参数
E = 3.45 * 10^3; % 弹性模量 (Pa)
I = 289.81; % 截面惯性矩 (m^4)
rho = 2500; % 梁材料密度 (kg/m^3)
A = 34.36; % 梁的截面积 (m^2)
L = 80; % 每跨长度 (m)
k = 2.96 * 10^9; % 弹簧支座刚度 (N/m)
n = 2; % 跨数
% 流体参数
Cm = 1.2; % 质量系数
rho_w = 1000; % 水的密度 (kg/m^3)
% 等效质量
m = rho * A; % 梁的单位长度质量
m_eff = m + (Cm - 1) * rho_w * A;
% 定义 alpha 作为 omega 的函数
alpha = (omega^2 * m_eff / (E * I))^(1/4);
% 定义符号矩阵A及分块矩阵
A_matrix = sym(zeros(4 * n, 4 * n)); % 确保A_matrix是符号矩阵
A1 = [0, 1, 0, 1;
0, -1, 0, 1];
A2 = [sin(alpha * L), cos(alpha * L), sinh(alpha * L), cosh(alpha * L);
cos(alpha * L), -sin(alpha * L), cosh(alpha * L), sinh(alpha * L);
-sin(alpha * L), -cos(alpha * L), sinh(alpha * L), cosh(alpha * L);
-cos(alpha * L) -(k / (alpha^3 * E * I)) * sin(alpha * L), sin(alpha * L) - (k / (alpha^3 * E * I)) * cos(alpha * L), cosh(alpha * L) - (k / (alpha^3 * E * I)) * sinh(alpha * L), sinh(alpha * L) - (k / (alpha^3 * E * I)) * cosh(alpha * L)];
A3 = [0, -1, 0, -1;
-1, 0, -1, 0;
0, -1, 0, 1;
1, 0, -1, 0];
A4 = [sin(alpha * L), cos(alpha * L), sinh(alpha * L), cosh(alpha * L);
-sin(alpha * L), -cos(alpha * L), sinh(alpha * L), cosh(alpha * L)];
% 构建A_matrix
for i = 1:n-1
% 更新A_matrix的分块
A_matrix(1:2, 1:4) = A1;
A_matrix(4*i-1:4*i+2, 4*i-3:4*i) = A2;
A_matrix(4*i+1:4*i+4, 4*i-3:4*i) = A3;
A_matrix(4*n-1:4*n, 4*n-3:4*n) = A4;
end
% 定义行列式的函数，使用subs将符号矩阵A_matrix的omega替换为数值omega_val
determinant_function = @(omega_val) double(det(subs(A_matrix, omega, omega_val)));
% 定义初始猜测值
omega_initial_guesses = 0.00001;
% 设置 fsolve 的精度选项
options = optimset('TolFun', 1e-8, 'TolX', 1e-8, 'MaxIter', 500, 'MaxFunEvals', 500);
% 求解
omega_solution = fsolve(determinant_function, omega_initial_guesses, options);

2025年01月02日 08点01分 6