计算科学 - 带有矩阵初始条件的 ode45 - 吾爱随笔录

带有矩阵初始条件的 ode45

计算科学 matlab 颂耦合

2021-11-29 05:12:46

编辑：我们有一个耦合系统，每个系统有 10 个颂歌。耦合出现在最后一个等式中。我考虑过使用 10 x 2 的矩阵作为初始条件。

我也在这里关注了一个具有相同标题的类似问题，但我仍然得到相同的错误（'输入参数太多。'）

time = [0 5];
x1_0 = [0 0 0 0 0 0 0 0 0 0];
x2_0 = [0 0 0 0 0 0 0 0 0 0];
initial = [x1_0;x2_0];
x = NaN(length(initial),2*length(time));

[t,x] = ode45(@ode,time,initial);


function [dxdt] = ode(x)
N = 2;
dxdt = NaN(10,2);
A = 3.25e-3;
B = 22e-3;
a = 100;
b = 50;
C = 135;
C1 = C;
C2 = 0.8*C;
C3 = 0.25*C;
C4 = 0.25*C;
C12 = C;

for i = 1:N
    dxdt(1,i) = x(6,i);
    dxdt(6,i) = A*a*C1*sigm((x(3,i)-x(4,i)+x(5,i))) - 2*a*x(6,i) - 
x(1,i)*a^2;
    dxdt(2,i) = x(7,i);
    dxdt(7,i) = A*a*C2*sigm((x(3,i)-x(4,i)+x(5,i))) - 2*a*x(7,i) - 
x(2,i)*a^2;
    dxdt(3,i) = x(8,i);
    dxdt(8,i) = A*a*C3*sigm((x(1,i))) - 2*a*x(8,i) - x(3,i)*a^2;
    dxdt(4,i) = x(9,i);
    dxdt(9,i) = B*b*C4*sigm(x(2,i)) - 2*b*x(9,i) - x(4,i)*b^2;
    dxdt(5,i) = x(10,i);
    if i == 1
        j = 2;
    elseif i == 2
        j = 1;
    end
    dxdt(10,i) = A*a*C12*sigm(x(3,j)-x(4,j)+x(5,j)) - x(10,i) - x(5,i);
end
end

function X = sigm(u)
u0 = 6e-3;
e0 = 2.5;
r = 0.56e3;
X = 2*e0/(1+exp(r*(u0-u)));
end

如果我的错误是使用矩阵初始条件而不是向量，使用 1 x 20 向量，并相应地调整 ode 形式是不切实际的，我认为 - 解决初始条件的另一种更有效的方法是什么 -我给出的输入是不必要的，为什么？- 是否有任何其他计算方式来表示耦合？

2个回答

您需要编写问题，使未知数是单个向量，而不是矩阵。在你的例子中 $N=2$ ，你将有一个未知的向量 $x(t)$ 大小的 $20\times 1$ （不是矩阵 $10\times 2$ ）。您将解决以下形状的问题

\dot{x} (t) = A (t) x (t), w i t h x \in R^{n}, A (t) \in R^{n \times n}

$\dot{x}(t) = A(t)\, x(t), \mathrm{ with }\ x \in \mathbb{R}^{n},\ A(t) \in \mathbb{R}^{n\times n}$ 在哪里

n = 190 \times 10

$n = 190 \times 10$ 在您的特定情况下。

如果您的示例确实代表了您的问题，则矩阵 $A(t)$ 将有不断的条目，所以 $A(t) \equiv A(0) \equiv A$ . 因此可以预先构建它（在集成循环之外）。在定义右侧（ode(x)在您的情况下）的函数中，您将只有一个矩阵向量乘法 $A\, x(t)$ . 并给出你的例子，你的矩阵 $A$ 将是稀疏的，因此您可以使用稀疏矩阵向量乘法从中受益。

PS。请注意，我认为 ode45 需要一个带有签名 (t,x) 的函数作为右侧。因此，您可能必须（正式）将您的函数编码为

function dxdt = ode(t,x)
    dxdt = A*x
end

下面的工作（从语法的角度，而不是模拟的角度）代码。我刚刚线性化了你的矩阵 $m\times n$ 带索引 $(i,j)$ 作为带有索引的向量 $(i+(j-1)*n)$ ; 重新排序功能，使用正确的签名ode45。

time = [0 5];
%x1_0 = [0 0 0 0 0 0 0 0 0 0];
%x2_0 = [0 0 0 0 0 0 0 0 0 0];
N = 2;
initial = zeros(10*N,1);
%x = NaN(length(initial),2*length(time));

function X = sigm(u)
u0 = 6e-3;
e0 = 2.5;
r = 0.56e3;
X = 2*e0/(1+exp(r*(u0-u)));
end



function dxdt = rhs(t,x)
  N=2;
  dxdt = NaN(10*N,1);
  A = 3.25e-3;
  B = 22e-3;
  a = 100;
  b = 50;
  C = 135;
  C1 = C;
  C2 = 0.8*C;
  C3 = 0.25*C;
  C4 = 0.25*C;
  C12 = C;

  for i = 1:N
      dxdt(1+(i-1)*10) = x(6+(i-1)*10);
      dxdt(6+(i-1)*10) = A*a*C1*sigm((x(3+(i-1)*10)-x(4+(i-1)*10)+x(5+(i-1)*10))) - 2*a*x(6+(i-1)*10) - x(1+(i-1)*10)*a^2;
      dxdt(2+(i-1)*10) = x(7+(i-1)*10);
      dxdt(7+(i-1)*10) = A*a*C2*sigm((x(3+(i-1)*10)-x(4+(i-1)*10)+x(5+(i-1)*10))) - 2*a*x(7+(i-1)*10) - x(2+(i-1)*10)*a^2;
      dxdt(3+(i-1)*10) = x(8+(i-1)*10);
      dxdt(8+(i-1)*10) = A*a*C3*sigm((x(1+(i-1)*10))) - 2*a*x(8+(i-1)*10) - x(3+(i-1)*10)*a^2;
      dxdt(4+(i-1)*10) = x(9+(i-1)*10);
      dxdt(9+(i-1)*10) = B*b*C4*sigm(x(2+(i-1)*10)) - 2*b*x(9+(i-1)*10) - x(4+(i-1)*10)*b^2;
      dxdt(5+(i-1)*10) = x(10+(i-1)*10);
      if i == 1
          j = 2;
      elseif i == 2
          j = 1;
      end
      dxdt(10+(i-1)*10) = A*a*C12*sigm(x(3+(j-1)*10)-x(4+(j-1)*10)+x(5+(j-1)*10)) - x(10+(i-1)*10) - x(5+(i-1)*10);
  end
end

[t,x] = ode45(@rhs,time,initial);

如果这不需要非常快，只需使用“reshape”函数在状态的矩阵和向量表示之间进行更改。这将导致更多的内存分配，这会大大减慢速度。

其它你可能感兴趣的问题

上一篇加权 QR 实现下一篇参考请求：Riks 方法（非线性 FEM）