Как получить больше выходов от ODE45 - PullRequest
Новая
       38

Как получить больше выходов от ODE45

0 голосов
/ 09 ноября 2019

У меня N = 2 связанных нелинейных динамических системы, связь которых задается матрицей 2 на 2 W. Каждая из них описывается n = 8 одой первого порядка. Приведенный ниже код решает эту связанную систему для многих значений параметра p: 

for i=1:length(p)
    [t,y(:,:,i)] = ode45(@(t,y) ode(t,y,n,N,W,p(i,:)), t, y0);
end

function [dydt] = ode(t,y,n,N,W,p)
    dydt = NaN(n, N);
    y = reshape(y,[n, N]);    
    y_out = zeros(N,1);
    F_Global = zeros(N,1);
    for i = 1:N
        y_out(i) = y(3,i)-y(4,i);
    end
    for i = 1:N
        F_Global(i) = W(i,:)*sigm(y_out);
    end
    for i = 1:N
        dydt(1,i) = y(5,i);
        dydt(2,i) = y(6,i);
        dydt(3,i) = y(7,i);
        dydt(4,i) = y(8,i);
        dydt(5,i) = sigm(y(3,i) - y(4,i)) - y(5,i) - y(1,i) + F_Global(i);
        dydt(6,i) = sigm(y(3,i) - y(4,i)) - y(6,i) - y(2,i);
        dydt(7,i) = C2*sigm(y(1,i)) + p(i) - y(7,i) - y(3,i);
        dydt(8,i) = sigm(y(2,i)) - y(8,i) - y(4,i);
    end
    dydt = reshape(dydt,n*N, 1);
end

function X = sigm(u)
    ...
end

Внутри функции уже вычислена разница: 

y_out(i) = y(3,i)-y(4,i);

Для i = 1, ..., N, и для всех времен и всех значений p предполагается, что это трехмерная матрица измерений 

y_out = size(length(time), length(p), N);

Кроме того, внутри функции вычисляется: 

F_Global(i) = W(i,:)*sigm(y_out);

, что для i = 1, ..., N и для всех значений p , но после усреднения по времени предполагается, что это двумерная матрицаразмеры 

F_Global = size(length(p),N);

Мне нужна помощь для извлечения y_out и F_Global в качестве выходных данных ode45

1 Ответ

0 голосов
/ 10 ноября 2019

измененный код 

output_potential = NaN(length(time),N,length(p));
Sigm = NaN(N,length(p));
input_firingRate = NaN(N,length(p));

for i=1:length(p)
    [t,y(:,:,i)] = ode45(@(t,y) ode(t,y,n,N,W,p(:,i)), time, y0);
    for j = 1:N
        output_potential(:,j,i) = y(:,3+(j-1)*n,i) - y(:,4+(j-1)*n,i);
        Sigm(j,i) = sigm(mean(output_potential(:,j,i)));
    end
    input_firingRate(:,i) = p(:,i) + A*a*W(:,:)*Sigm(:,i);
end
...