1

我是 MATLAB 的新手用户。我有代码试图找到状态空间模型的时间历史。有四个一阶 ODE，我想同时使用ode45. 待解方程的实质如下：

x1_dot = x2

x2_dot = -[M] * [K] * x1 - [M] * [C] * x2 + constant*[M] * [P3] * x3 + constant*[M] * [P4] * x4

x3_dot = x2 - constant*x3

x4_dot = x2 - constant*x4

其中[M], [K], [C], [P3], 和[P4]是 3x3 矩阵；x1, x2, x3,x4都是 3x1 向量；等x1_dot表示时间导数（它们是 3x1 向量）。我只有初始条件x1。

我编写的 MATLAB 代码如下。这段代码在我的整个程序中。我没有调用单独的函数，因为我不知道如何通过ode45函数传递所有矩阵/向量。我收到错误消息：“索引超出矩阵维度。”

tspan = 0:1:20;

initial = [0 0.03491 0];

f = @(t,x) [x(2);
            -inv(M_Dbl_Bar_Matrix)*K_Dbl_Bar_Matrix*x(1) - inv(M_Dbl_Bar_Matrix)*C_Dbl_Bar_Matrix*x(2) + (0.5*rho*U^2)*inv(M_Dbl_Bar_Matrix)*P3_Matrix*x(3) + (0.5*rho*U^2)*inv(M_Dbl_Bar_Matrix)*P4_Matrix*x(4);
            x(2) - Beta_1*x(3);
            x(2) - Beta_2*x(4)];

[t,xp] = ode45(f,tspan,initial);

问题：

  1. 如何处理x(1)、x(2)、x(3)和x(4)inode45是 3x1 向量？

  2. 如何为这个方程组应用初始条件？例如，我是否使用矢量格式，例如： initial = [0 0.03491 0; 0 0 0; 0 0 0; 0 0 0]？

  3. 我是否正确使用/编写了函数（f）和ode45？

4

2 回答 2

0

只需将 12×1 系统集成为 12 个耦合 ODE。

其他一些观察：

  1. inv()尽可能避免- 它缓慢且不准确。请改用mldivide（或mrdivide）。此外，您在 的每次评估中重新计算不少于4次f，而它基本上是恒定的！
  2. 您似乎希望每秒都有输出。ode45是一个可变步长积分器，这意味着它会根据对每一步误差的估计自动调整步长。在偏离所选择的时间的时间请求输出ode45（称为“密集输出”）很容易做到，但它会花费你额外的函数评估。通常这并不是真正需要的，您可以改用更有效的方法tspan = [0 20]。这完全取决于您的具体需求。

现在，这就是我想出的：

% Time interval of interest
tspan = [0 20];

% Initial values
x1_0 = [0 0.03491 0];
x2_0 = [0       0 0];
x3_0 = [0       0 0];
x4_0 = [0       0 0];
x0   = [x1_0 x2_0 x3_0 x4_0].';

% Pre-compute a few constants
Fd  = 0.5*rho*U^2;
P3f = Fd*M_Dbl_Bar_Matrix\P3_Matrix;
P4f = Fd*M_Dbl_Bar_Matrix\P4_Matrix;
P1f = -M_Dbl_Bar_Matrix\K_Dbl_Bar_Matrix;
P2f = -M_Dbl_Bar_Matrix\C_Dbl_Bar_Matrix;

% The derivative (12×1, but constructed as 4·3×1)

one = 1:3;   three = 7:9;   % well-named index vectors to
two = 4:6;   four  = 10:12; % make our lives a bit easier

f = @(t,x) [x(two)
            P1f*x(one) + P2f*x(two) + P3f*x(three) + P4f*x(four)
            x(two) - Beta_1*x(three)
            x(two) - Beta_2*x(three)];

Now integrate this system
[t, xR] = ode45(f, tspan, x0);

% extract results in the same kind of blocks:
x1 = xR(:, one);
x2 = xR(:, two);
x3 = xR(:, three);
x4 = xR(:, four);

% ... process the results in whatever way you see fit
于 2017-01-05T09:45:52.830 回答
0

1.) “成为 3x1 向量”是什么意思？x(i)必须是单个变量才能使其工作，因此 size(x)=1x4。在 ODE 的上下文中，拥有x(i)=(x,y,z)并没有任何意义。

2.) 您的变量向量x的长度应为 4，因此您的初始条件应反映这一点。

initial = [0 0.03491 0 0];

对我来说很好。

3.) 是的，我想是的。

您能否提供有关您正在尝试做的事情的更多信息？

编辑

好的，我想我现在明白你想要做什么了。

例如，您有 4 个向量：x(x,y,z)、x'(x,y,z)、p(x,y,z)、q(x,y,z) 和一个大矩阵 4x4 矩阵3x3 矩阵。所以我猜

（原谅托管的图像，因为我是一个 stackoverflow 菜鸟，我不允许直接发布 html 或图像，或者多个链接）

数学 1、2、3 参考：https ://postimg.org/gallery/1qh2ywiqq/

来自链接的数学 1

和

来自链接的数学 2

你的状态空间方程。

所以要解决这个问题，你必须设置一个 12 维 ode45 问题，而不是现在的 4 维问题，因为每个向量都有 3 个分量。您需要做的是通过明确指定每个条目将 4x4 矩阵更改为 12x12 矩阵。您还必须给出f = @(t,x)一个 1x12 向量：

来自链接的数学 3

然后设置initial为初始条件的 1x12 向量（向量中的每个维度一个）。

有了我们现在拥有的矩阵形式和初始条件，我们可以使用这个来源：https ://nl.mathworks.com/help/symbolic/solve-a-system-of-differential-equations.html#buxuujb

它告诉我们如何正确设置它。

我目前没有时间给你代码，但我希望我能正确理解你想要做什么，这将是有用的。

于 2017-01-04T21:56:32.050 回答