我正在尝试在 MATLAB 中求解达芬方程。
在哪里
为此，我编写了一个用于 ode45 的函数。
function [F] = ode_fun(t,y)
F = zeros(size(y));
F(1) = y(2);
F(2) = -k*y(1)-c*y(2)-k3*y(1)^3+A*sin(omega*t);
end
然后我调用 ode 函数：
[t,u] = ode45('ode_fun',[0 T],[0 0]);
这是正确的方法吗？
因为我希望在绘制“力与位移”图时获得三次函数，但实际上我获得了一条直线。
编辑：我也会给你我用于模拟的数值。
基本上，我从一个示例中找到的线性方程开始，然后添加了非线性项来检查我的代码是否能够区分线性和非线性动态系统。
线性情况的值是：,,.
然后我将非线性项添加为
作为强制项，我决定使用正弦波和.
最终模拟时间为
编辑：我正在测试的完整代码是我在另一个问题中谈到的，我从 3D 表面发布了 2D 横截面
我试图了解此过程是否适用于非线性参数识别
在某些情况下，三次方看起来非常线性。您的线性刚度参数，和, 相当高。如果您查看具有指定值的位置状态变量的范围，您会注意到系统仅偏离零大约. 从物理上讲，这实际上是在系统的线性状态中，因此任何非线性分量都将非常小。减少通过一个数量级，振荡更大，并且在系统进入稳态后，力（加速度）与位置的关系图开始看起来是立方的：
一些建议。
1）我可能是一个好主意，开始为这个系统使用一个刚性求解器，这取决于你的参数是什么。在Matlab中，ode15s是一个不错的选择。
2）在你的另一个问题中查看你的代码，你真的应该将你的参数值作为变量传递，而不是在多个地方硬编码它们。有关详细信息，请参阅这篇关于Matlab 中的参数化函数的文章。对于您的代码：
tspan = [0 60];
k = 1e3;
c = 40;
k3 = 1.1e5;
A = 30;
omega = 10;
ode_fun = @(t,u,k,c,k3,A,omega)[u(2);-k*u(1)-c*u(2)-k3*u(1).^3+A*sin(omega*t)];
[t,u] = ode15s(@(t,u)ode_fun(t,u,k,c,k3,A,omega),tspan,[0 0]);
figure;
plot(t,u);
xlabel('Time');
ylabel('State');
acc = -k*u(:,1)-c*u(:,2)-k3*u(:,1).^3+A*sin(omega*t);
figure;
plot(u(:,1),acc);
xlabel('Position');
ylabel('Acceleration');
在这里，我使用了 ODE 函数直接获取加速度。
3) 您可能希望将积分时间 ( tspan) 指定为具有两个以上元素的向量，以确保特定的输出步长。（注意，这不是积分器内部使用的实际步长 -请参见此处。）如果您使用有限差分方案来计算加速度，这将有所帮助，并且在某些情况下它可以产生更平滑的图。