Filter löschen
Filter löschen

Adjacency matrix for ODE45

8 Ansichten (letzte 30 Tage)
N/A
N/A am 4 Sep. 2024
Kommentiert: N/A am 5 Sep. 2024
Hi
I am trying to model a systema as below:
where kij is :
I tried to model two sets of states manually first, without Aij (the adjacency matrix) as below:
clear all;clc;
%%
global b r I x0 alpha beta gamma
b=3.0;
r=0.006;
I=2.8;
x0=-1.6;
alpha=1;
beta=10;
gamma=0.5;
initial_condition=[-1,0.2,0.4,0.8,0,0.2,0.1,0.2];
tspan=[0 3000];
[t,y]=ode45(@EquationSys,tspan,initial_condition);
plot(t,y(:,1))
hold on
%%
function dy=EquationSys(t,y)
global b r I x0 alpha beta gamma
x1=y(1); %-1
y1=y(2);%0.2
z1=y(3);%0.4
x2=y(4);%0.8
y2=y(5);%0
z2=y(6);%0.2
k12=y(7);%0.1
k21=y(8);%0.2
dy=[y1-x1^3+b*x1^2-z1+I+k12*(x2-x1);
1-5*x1^2-y1;
r*(4*(x1-x0)-z1);
y2-x2^3+b*x2^2-z2+I+k21*(x1-x2);
1-5*x2^2-y2;
r*(4*(x2-x0)-z2);
k12*(alpha*exp(-beta*(x1 - x2)^2) - gamma*(k12+1));
k21*(alpha*exp(-beta*(x2 - x1)^2) - gamma*(k21+1));
];
end
%
%%%%%%%%%%%%%%%%%
For some reason, I want to add Aij to this model and vectorization the code and I programmed the code as below:
clear all;clc;
% HR Parameters
global b r x0 I N alpha beta gamma Adj_Matrix
b=3.0;
r=0.006;
I=2.8;
x0=-1.6;
alpha=1;
beta=10;
gamma=0.5;
%%
% Simulation parameters
T = 3000; % Total time
time_span = [0 T]; % Time vector
%%
% Number of neurons
% n=2; %Number of neurons on each edges
N = 2; % Number of all neurons
%% Adjacency Matrix
% Initialize the matrix with zeros
Adj_Matrix = zeros(N, N);
Adj_Matrix(1,2)=1;
Adj_Matrix(2,1)=1;
%%
% Flatten initial conditions and initial coupling strengths into a single vector
initial_conditions = [-1,0.8,0.2,0,0.4,0.2,0,0.1,0.2,0];
%%
% Solve the ODE for the coupled system with STDP
[time, sol] = ode45(@ hr_equations, time_span, initial_conditions);
%x, y, z, and k from the solution
x = sol(:, 1:N);
y = sol(:, N+1:2*N);
z = sol(:, 2*N+1:3*N);
k = sol(:, 3*N+1:end);
%%
% % Plotting results
%
% % If you want to visualize all the neurons' behavior together
plot(time,x(:,1))
hold on
%%
function dstate_dt = hr_equations(time, state)
global b r x0 I N alpha beta gamma Adj_Matrix
% Reshape the state vector
x = state(1:N);
y = state(N+1:2*N);
z = state(2*N+1:3*N);
k = reshape(state(3*N+1:end), N, N);
% HR Model Equations
dx = y-x.^3+b*x.^2-z+I+sum((Adj_Matrix.*k).*(x'-x),2);
dy = 1-5*x.^2-y;
dz = r*(4*(x-x0)-z);
dk = zeros(N, N);
for i = 1:N
for j = 1:N
if i ~= j
dk(i,j) = k(i,j) * (alpha * exp(-beta * (x(i) - x(j))^2) - gamma * (k(i,j) + 1));
end
end
end
% Flatten the derivatives to return as a single column vector
dstate_dt = [dx; dy; dz; reshape(dk, [], 1)];
end
However, after running this code and comparing two plots I found that there is slight differences between two results. I cannot find the issue. Can you help me in this regard?
Thanks

Akzeptierte Antwort

Torsten
Torsten am 4 Sep. 2024
Verschoben: Torsten am 4 Sep. 2024
k = [0,0.1,0.2,0];
reshape(k,2,2)
ans = 2x2
0 0.2000 0.1000 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
I think you want the reshaped matrix to be
reshape(k,2,2).'
ans = 2x2
0 0.1000 0.2000 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

Weitere Antworten (0)

Kategorien

Mehr zu Particle & Nuclear Physics finden Sie in Help Center und File Exchange

Produkte


Version

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by