如何在Matlab中求解耦合系统的初始条件

Question

我们有一个10个ode的耦合系统。耦合出现在最后一个方程中。

我考虑使用一个10乘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);

for i = 1:N
    dxdt(1,i) = x(6,i);
    dxdt(6,i) = (x(3,i)-x(4,i)+x(5,i)) - x(6,i) - x(1,i);
    dxdt(2,i) = x(7,i);
    dxdt(7,i) = (x(3,i)-x(4,i)+x(5,i)) - x(7,i) - x(2,i);
    dxdt(3,i) = x(8,i);
    dxdt(8,i) = (x(1,i)) - x(8,i) - x(3,i);
    dxdt(4,i) = x(9,i);
    dxdt(9,i) = (x(2,i)) - x(9,i) - x(4,i);
    dxdt(5,i) = x(10,i);
    if i == 1
        j = 2;
    elseif i == 2
        j = 1 ;
    end
    dxdt(10,i) = (x(3,j)-x(4,j)+x(5,j)) - x(10,i) - x(5,i);

end
end

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

EDIT:如果我使用选项："function dxdt = ode(t，x)“，错误如下：

Error in Test11>ode (line 33)
dxdt(10,i) = (x(3,j)-x(4,j)+x(5,j)) - x(10,i) - x(5,i);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, 
varargin);

Error in Test11 (line 9)
[t,x] = ode45(@ode,time,initial);

如果我改变调用函数的方式，那么："t，x= ode45(@(t，x)ode(x)，time，initial);“，结果如下：

Not enough input arguments.

Error in Test11>ode (line 19)
dxdt(1,i) = x(6,i);

Error in Test11>@(t,x)ode(x)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, 
varargin);

Error in Test11 (line 9)
[t,x] = ode45(@(t,x)ode(x), time, initial);

Answer 1

错误'Too many input arguments.'是因为ode函数只接受一个变量x，而ode45向它传递了两个变量，时间实例t和状态x。即使您的常微分方程与当前时间无关，您也需要处理此输入参数。

让您的ode函数接受两个输入参数：

function dxdt = ode(t,x)
% etc

或者，更改调用ode函数的方式：

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

此外，我不认为ode45可以处理“状态矩阵”，所以我确实会将垂直于一个大状态向量的两个状态向量连接起来。如果你不想重写你的代码，在循环前后添加一些重塑：

function [dxdt] = ode(xv)
    N = 2;
    dxdt = NaN(10,2);

    % reshape to two column vectors
    x = reshape(xv,[],2); % xv is one large column vector

    for ki = 1:N
        dxdt(1,ki) = x(6,ki);
        dxdt(6,ki) = (x(3,ki)-x(4,ki)+x(5,ki)) - x(6,ki) - x(1,ki);
        dxdt(2,ki) = x(7,ki);
        dxdt(7,ki) = (x(3,ki)-x(4,ki)+x(5,ki)) - x(7,ki) - x(2,ki);
        dxdt(3,ki) = x(8,ki);
        dxdt(8,ki) = (x(1,ki)) - x(8,ki) - x(3,ki);
        dxdt(4,ki) = x(9,ki);
        dxdt(9,ki) = (x(2,ki)) - x(9,ki) - x(4,ki);
        dxdt(5,ki) = x(10,ki);
        if ki == 1
            kj = 2;
        elseif ki == 2
            kj = 1;
        end
        dxdt(10,ki) = (x(3,kj)-x(4,kj)+x(5,kj)) - x(10,ki) - x(5,ki);
    end

    % reshape dxdt to one column vector
    dxdt = reshape(dxdt,[],1);

end