带有矩阵初始条件的 ode45答案

【问题标题】：ode45 with matrix initial conditions带有矩阵初始条件的 ode45
【发布时间】：2017-05-19 06:59:17
【问题描述】：

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

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);

问题：

如何处理 ode45 中的 x(1)、x(2)、x(3) 和 x(4) 作为 3x1 向量？
如何为这个方程组应用初始条件？例如，我是否使用矢量格式，例如：initial = [0 0.03491 0; 0 0 0; 0 0 0; 0 0 0]？
我是否正确使用/编写了函数 (f) 和 ode45？

【问题讨论】：

标签： matlab ode45

【解决方案1】：

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) 和一个大矩阵由 3x3 矩阵组成的 4x4 矩阵。所以我猜

（请原谅托管的图片，因为我是 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

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

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

【讨论】：

你遇到了这个问题。我正在求解一个矩阵系统：{x-dot} = [A]{x}。 {x} 由四个 3x1 向量（x1、x2、x3、x4）组成，产生一个整体 12x1 向量。 [A] 是一个由四行四列组成的矩阵，其中每个行列条目是一个 3x3 矩阵；所以 [A] 总共是一个 12x12 的矩阵。这个矩阵方程可以写成我上面提到的四个一阶 ODE。每个 {x} 向量都有初始条件，所以我应该有 initial = transpose([0 0.03491 0 0 0 0 0 0 0 0 0 0])。这是一个 12x1 的初始条件向量。这个问题应该是ode45解决的，不知道怎么解决的。

【解决方案2】：

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

其他一些观察：

尽可能避免使用inv()——它速度慢且不准确。请改用mldivide（或mrdivide）。此外，在每次评估 f 时，您都会重新计算不少于 4 次，而它基本上是恒定的！
您似乎希望每秒都有输出。 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

【讨论】：