DQM法求解微分方程

求解方程 𝑦′′(𝜓,𝑡)+𝑎1𝑦′(𝜓,𝑡)+𝑎2𝑦(𝜓,𝑡)+𝑎3cos⁡(𝑎4𝑡)𝑦(𝜓,𝑡)=0y′′(ψ,t)+a1​y′(ψ,t)+a2​y(ψ,t)+a3​cos(a4​t)y(ψ,t)=0 在 𝜓=0ψ=0 处的动力系统响应

1. DQM矩阵的构建：在构建一阶和二阶导数近似矩阵时，如果处理不当，可能会导致矩阵尺寸不匹配。
2. 更新解的循环：在更新解的循环中，如果对矩阵的操作不正确，也可能导致尺寸不匹配。

以下是运行的代码，特别注意矩阵尺寸的匹配

% 定义常数 a1 = 1;

% 线性阻尼系数 a2 = 2;

% 二次刚度系数 a3 = 3;

% 非线性系数 a4 = 4; % 频率

% 初始条件

y0 = 0.00001;

% y(0) = 0.00001

yp0 = 0.000001;

% yp(0) = 0.000001

% 网格数量 N = 20;

% 网格点数量

% 切比雪夫节点

x = cos((0:N-1)*(2*pi/(N-1)));

% 初始化微分求积近似矩阵

D0 = eye(N);

% 零阶导数近似矩阵

D1 = zeros(N, N);

% 一阶导数近似矩阵

D2 = zeros(N, N);

% 二阶导数近似矩阵

% 构建一阶导数近似矩阵 D1

for i = 1:N

for j = 1:N

if i ~= j D1(i, j) = -1 / (x(i) - x(j));

else D1(i, j) = 0;

end

end

end

% 构建二阶导数近似矩阵 D2

D2 = D1 * D1';

% 应用边界条件

y(-1, t) = y(1, t) = 0

D1(1, :) = 0;

D1(end, :) = 0;

D2(1, :) = 0;

D2(end, :) = 0;

% 时间跨度

T = 10;

% 总时间

dt = 0.01;

% 时间步长

t_steps = round(T / dt);

% 初始化解矩阵

Y = zeros(N, t_steps+1);

Y(:, 1) = y0 * ones(N, 1);

% 初始化y(psi,0) % 时间迭代求解

for t = 0:t_steps-1 t_actual = (t + 1) * dt;

% 非齐次项 nonHomogeneousTerm = a3 * cos(a4 * t_actual);

% 更新y的值

Y(:, t+2) = Y(:, t+1) + dt^2 * (-a1 * D1 * Y(:, t+1) - a2 * D2 * Y(:, t+1) - nonHomogeneousTerm * D0 * Y(:, t+1));

end

% 绘制结果

time = (0:t_steps) * dt; plot(time, Y(round(end/2)+1, :), 'b-');

% 绘制中间节点的响应

xlabel('Time (t)');

ylabel('Displacement (y)');

title('Dynamic Response of the System at Psi = 0 using DQM');

我们首先定义了问题的常数和初始条件，然后创建了切比雪夫节点 x 和微分求积近似矩阵 D0、D1 和 D2。我们手动构建了一阶导数近似矩阵 D1 和二阶导数近似矩阵 D2。

接下来，我们应用了边界条件，并初始化了解矩阵 Y。在时间迭代求解中，我们更新了 Y。

最后，我们绘制了中间节点（即 𝜓=0ψ=0 处）的系统响应随时间的变化。